EMENTAL  ANALYSIS  OF  NONLINEAR  STRUCTURAL  MECHANICS  PROBLEMS  WITH  APPLICATIONS  TO  GENERAL 

SHELL  STRUCTURES  R^m*_4J42 


gvtdoc 

D  211, 
9: 

4142 


1  SHIP  RESEARCH  AND  DEVELOPMENT  CENTER 

)  Bethesda,  Maryland  20034 


NIMITZ  LIBRARY  -  U.S.  NAVAL  ACADEMY 

3 1  9337  00373 ' 7005 


INCREMENTAL  ANALYSIS  OF  NONLINEAR  STRUCTURAL 
MECHANICS  PROBLEMS  WITH  APPLICATIONS  TO 
GENERAL  SHELL  STRUCTURES 


by 

Rembert  F.  Jones,  Jr. 


APPROVED  FOR  PUBLIC  RELEASE: 
DISTRIBUTION  UNLIMITED 


STRUCTURES  DEPARTMENT 
RESEARCH  AND  DEVELOPMENT  REPORT 


September  1973 


DEC  5 


Report  4142 


'The  Nova!  Ship  Research  and.  Development  Center  is  a  U.  S.  Navy  center  for  laboratory 
effort  directed  at  achieving  improved  sea  and  air  vehicles.  It  was  formed  in  March  1967  by 
merging  the  David  Taylor  Model  Basin  at  Carderocfc,  Maryland  with  the  Marine  Engineering 
Laboratory  at  Annapolis,  Maryland* 


Naval  Ship  Research  and  Development  Center 
Betheada,  Md.  20034 


MAJOR  NSRDC  ORGANIZATIONAL  COMPONENTS 


NDW-NSRDC  3960/44  (REV.  8/71) 

GPO  917-872 


DEPARTMENT  OF  THE  NAVY 
NAVAL  SHIP  RESEARCH  AND  DEVELOPMENT  CENTER 

BETHESDA,  MARYLAND  20034 


INCREMENTAL  ANALYSIS  OF  NONLINEAR  STRUCTURAL 

MECHANICS  PROBLEMS  WITH  APPLICATIONS  TO 
GENERAL  SHELL  STRUCTURES 


by 

Rembert  F.  Jobes,  Jr. 


APPROVED  FOR  PUBLIC  RELEASE: 
DISTRIBUTION  UNLIMITED 


September  1973 


Report  4142 


TABLE  OF  CONTENTS 


Page 


ABSTRACT .  1 

ADMINISTRATIVE  INFORMATION  .  1 


Chapter 

1.  INTRODUCTION . 

1.1  GENERAL  REMARKS . 

1.2  NONLINEAR  ANALYSES  . . . 

1.2.1  Nonlinear  Material  Behavior . 

1.2.2  Large  Deformations . 

1.3  GENERAL  SHELL  FINITE  ELEMENTS . 

2.  GENERAL  INCREMENTAL  FORMULATION . 

2.1  DEFINITION  OF  STRAIN  AND  STRESS . 

2.2  PRINCIPLE  OF  VIRTUAL  WORK  -  INCREMENTAL  THEORY 

2.3  FINITE  ELEMENT  FORMULATION  OF  THE  EQUILIBRIUM 

EQUATION  . 

2.4  INITIAL  STRESS  PROBLEMS  . 

3.  CONSTITUTIVE  LAWS  OF  PLASTICITY . 

3.1  BASIC  ASSUMPTIONS  -  INITIAL  YIELD  CONDITION  . 

3.2  LOADING  CRITERIA  . 

3.3  WORK  HARDENING  -  THE  DRUCKER  POSTULATE  . 

3.4  FLOW  RULE  FOR  WORK-HARDENING  MATERIALS  . 

3.5  SUBSEQUENT  LOADING  SURFACES  -  HARDENING  RULES. 

3.6  EQUIVALENT  STRESS  AND  EQUIVALENT  STRAIN . 

3.7  STRESS-TO-STRAIN  TRANSFORMATION . 

3.8  STRAIN -TO-STRESS  TRANSFORMATION . 

4.  GEOMETRIC  DETAILS  . 

4.1  SURFACE  DEFINITION . 

4.2  ELEMENT  DEFINITION  AND  BOUNDARY . 

4.3  SLOPE  OF  SURFACE  IN  ELEMENT  SYSTEM . 

4.4  SURFACE  VECTORS  . 

5.  QUADRILATERAL  SHELL  FINITE  ELEMENT  . 

5.1  STRAIN  DISPLACEMENT  RELATIONSHIP  . 

5.2  ASSUMED  DISPLACEMENTS  . i . 

5.3  GENERALIZED  DISPLACEMENTS  . . . 

5.4  GENERALIZED  STRESSES  AND  STRAINS . 

5.5  STIFFNESS  MATRIX  AND  EQUILIBRIUM  CORRECTION 

VECTOR  . . 

5.6  LOADING  . 

5.7  TRANSFORMATIONS  . . . 


2 

2 

2 

2 

3 

5 

9 

9 

10 

16 

20 

21 

21 

22 

23 

23 

25 

25 

27 

30 

33 

33 

33 

35 

37 

40 

40 

43 

46 

50 

51 
51 
56 


» 


* 


n 


TABLE  OF  CONTENTS  (Contd.) 

Page 

5.8  CONVERGENCE  REQUIREMENTS .  60 

6.  IMPLEMENTATION  OF  THEORY .  63 

6.1  LINEAR  ELASTIC  LARGE  DEFORMATION  SOLUTIONS  .  63 

6.2  INELASTIC  LARGE  DEFORMATION  SOLUTIONS .  66 

7.  NUMERICAL  EXAMPLES  .  72 

7.1  RIGID  BODY  MODES .  72 

7.2  MEMBRANE  CYLINDER  -  STRESSES .  72 

7.3  GRAVITY-LOADED  CYLINDER  .  74 

7.4  RING-STIFFENED  CYLINDER .  77 

7.5  BAR-SPRING  PROBLEM  .  79 

7.6  SIMPLY  SUPPORTED  SQUARE  PLATE .  81 

7.7  CENTRALLY  LOADED  SPHERICAL  CAP .  83 

7.8  PEAR-SHAPED  CYLINDER .  83 

8.  SUMMARY  AND  RECOMMENDATIONS  FOR  FURTHER  WORK .  88 

8.1  SUMMARY .  88 

8.2  RECOMMENDATIONS  FOR  FURTHER  WORK  .  89 

ACKNOWLEDGMENTS  .  92 

APPENDIX  -  QUADRATIC  SURFACES  .  93 

REFERENCES . 94 


LIST  OF  FIGURES 


Figure 

2.1  Original  and  Deformed  Configurations  . 

3.1  Isotropic  Hardening  . 

3.2  Kinematic  Hardening  . . . 

3.3  Idealization  of  Plastic  Stress-Strain  Curve  . 

4.1  Coordinate  Systems . 

4.2  Element  in  Local  System  . 

4.3  Surface  Vectors  . 

5.1  Mapping  of  Shell  Surface  . 

5.2  Auxiliary  Coordinate  Systems . 

5.3  Differential  Area  on  Shell  Surface  and  in  Base  Plane . 

5.4  Element  Interface  on  Shell  Surface . 

6. 1  Region  of  Integration  . 

6.2  Section  Through  Element  . 

6.3  Calculation  of  Ratio  m  and  Normal  to  Yield  Surface  . 

7.1  Eigenvalues  of  Stiffness  Matrix  of  Quadrilateral  Element 

for  Sector  of  Sphere . 

7.2  Membrane  Cylinder  . 

7.3  Cylindrical  Shell  Roof . 

7.4  Total  Number  of  Degrees  of  Freedom  (Unknowns)  . 

7.5  10-Element  Idealization  of  Ring-Stiffened  Cylinder . 

7.6  One-Dimensional  Bar-Spring  Problem . 

7.7  Load  Deflection  Curve  for  Square  Plate  . 

7.8  Centrally  Loaded  Spherical  Cap  . 

7.9  Comparison  of  Computed  and  Experimental  Results  for 

Centrally  Loaded  Spherical  Cap  . 

7.10  Load -Displacement  for  Pear-Shaped  Cylinder . 

8.1  Nonlinear  Geometric  Pseudo-Force  Solution  . 

LIST  OF  TABLES 

Table 

5.1  B1  and  B2  Matrices  . 

7.1  Bar-Spring  Problem  -  Displacement  Versus  Tolerance  Ratio 


Page  « 

11 

26  , 

26 
31 
34 
34 
38 
41 
45 
53 
62 
64 
68 
70 

73 

75 

75 

76 
78 
80 
82 

84 

85 

86 
91 


52 

81 


» 


iv 


A 


ABSTRACT 


The  primary  theme  is  the  nonlinear  behavior  of  shell  structures.  A 
general  formulation  that  can  be  applied  to  all  structures  is  derived  and  then 
specialized  for  application  to  shells  through  the  development  of  a  quadrilateral 
shell  finite  element. 

The  formulation  is  in  incremental  form  and  is  derived  in  terms  of  the 
material  coordinate  system.  The  procedure  evolved  to  solve  the  nonlinear 
equilibrium  equations  allows  for  the  inclusion  of  all  terms  as  opposed  to  pre¬ 
vious  incremental  procedures  where  the  equations  are  linearized.  The 
derivation  of  the  formulation  and  solution  procedure  is  independent  of  the 
finite  element  method,  but  it  is  this  method  that  is  used  to  exercise  and  apply 
the  formulation.  An  elementary  nonlinear  problem  is  first  solved  and  then 
more  complex  problems  are  studied  by  using  the  formulation  and  solution 
procedure. 

The  quadrilateral  shell  element  which  is  introduced  uses  the  exact 
geometry  of  the  shell  surface  in  the  computation  of  stiffness  matrices  and 
loads.  A  consistent  method  is  used  for  transforming  the  element  and  its 
properties;  this  improves  element  effectiveness  for  application  to  branched  and 
reinforced  shells.  The  shell  surface  geometry  is  used  in  the  enforcement  of 
the  continuity  of  displacements  for  satisfying  the  convergence  requirements. 
The  ability  of  the  element  to  undergo  rigid  body  motion  without  causing 
strains  is  examined  by  studying  the  eigenvalues  of  the  stiffness  matrix  of  an 
element. 


ADMINISTRATIVE  INFORMATION 

The  study  was  conducted  at  the  Naval  Ship  Research  and  Development  Center  (NSRDC) 
under  the  sponsorship  of  the  Naval  Ships  Systems  Command.  Funding  was  provided  by 
Subproject  SF  43.422.210,  Task  15053,  Work  Unit  1-1725-395.  This  report  is  based  on  the 
dissertation  submitted  to  the  faculty  of  The  Catholic  University  of  America  in  partial  ful¬ 
fillment  of  the  requirements  for  the  Doctor  of  Philosophy  degree.  The  format  is  that  of 
Catholic  University  rather  than  NSRDC. 


1 
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CHAPTER  1 
INTRODUCTION 


1.1  GENERAL  REMARKS 

The  mathematical  description  of  large  finite  deformations  is  well  founded.  Unfor¬ 
tunately,  it  is  also  very  complicated  and  difficult  to  apply  to  complex  practical  problems. 
Writers  have  presented  the  elegant  theory  for  finite  deformations  but  have  been  hampered  in 
the  solution  of  other  than  very  basic  problems.  The  present  availability  of  high-speed  com¬ 
puters,  coupled  with  numerical  methods,  now  makes  it  feasible  to  exploit  the  significant 
fundamental  work  that  has  been  accomplished  by  pioneers  in  the  field. 

The  most  prominent  numerical  approach  for  solving  complex  structural  problems  is 
the  finite  element  method.  There  has  been  a  great  deal  of  development  and  application  of 
the  method  within  the  last  15  years.  More  investigators  in  the  solid  mechanics  and  structural 
engineering  fields  are  probably  devoted  to  this  subject  than  to  any  other.  Significant  progress 
has  been  made.  The  requirements  that  must  be  satisfied  to  achieve  adequate  solutions  have 
been  recognized  and  methods  have  been  conceived  to  satisfy  a  large  number  of  these  require¬ 
ments.  Progress  on  the  most  difficult  of  problems  has  been  significant. 

The  need  for  a  comprehensive  analysis  of  general  shell  structures  prompted  investi¬ 
gators  to  turn  to  the  finite  element  method.  The  early  work  was  directed  to  developing  linear 
analyses  for  shells,  which  is  now  well  in  hand.  Small  deflection,  elastic  analyses  of  shells 
with  arbitrary  geometry  and  reinforcement  can  now  be  conducted  routinely.  Research  is 
currently  being  conducted  to  extend  these  shell  analyses  into  the  nonlinear  geometry  and  non¬ 
linear  material  range.  Significant  work  has  already  been  done  but  reasonable  accuracy  is 
obtained  at  the  expense  of  very  high  computing  costs. 

The  research  reported  herein  is  directed  to  a  better  understanding  of  the  nonlinear 
behavior  of  structures,  in  particular,  that  of  shell  structures.  A  general  formulation  for  con¬ 
ducting  nonlinear  incremental  analyses  is  introduced  and  a  quadrilateral  shell  finite  element  is 
developed.  It  is  hoped  that  both  will  contribute  to  the  more  efficient  solution  of  nonlinear 
structural  problems. 

1.2  NONLINEAR  ANALYSES 

The  evolution  of  the  procedures  for  handling  material  and  geometrical  nonlinearities 
by  the  finite  element  method  has  followed  two  separate  and  distinct  paths.  This  has  been 
possible  because  the  formulation  of  one  type  of  nonlinearity  does  not  depend  on  the  other 
and  their  effects  can  be  included  independently. 

1.2.1  Nonlinear  Material  Behavior 

Two  successful  techniques  have  been  used  to  include  the  effects  of  time-independent 
isothermal,  small  strain  plasticity  behavior.  Both  are  based  on  the  incremental  theory  of 
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plasticity  which  has  more  physical  meaning  than  the  deformation  or  total  strain  theory.  The 
classical  approaches  to  solving  plasticity  problems  were  primarily  restricted  to  using  the  defor¬ 
mation  theory,  but  this  is  no  longer  necessary. 

The  first  approach  for  solving  plasticity  problems  numerically  takes  advantage  of  the 
analogy  between  plastic  strain  and  body  forces  first  recognized  by  Ilyushin  [1  ]*.  The  total 
strain  is  decomposed  into  the  elastic  and  the  plastic  portions.  Element  stiffness  matrices  are 
computed  on  the  assumption  of  full  elastic  behavior  in  the  element  while  the  plastic  strains 
are  converted  to  equivalent  pseudo  forces  at  the  nodes.  The  plastic  strain  increments  used  in 
the  plastic  load  vector  are  those  that  were  computed  in  the  previous  increment.  Applications 
of  this  approach  are  described  by  Armen  et  al.  [2]  and  Argyris  et  al.  [3). 

The  second  approach  used  to  incorporate  inelastic  material  behavior  is  the  tangent 
stiffness  method.  This  is  the  formulation  that  is  used  in  the  present  study  and  will  be  de¬ 
scribed  in  detail  subsequently.  Briefly,  the  element  stiffness  matrices  are  recomputed  at  each 
load  step  by  using  the  current  transformation  of  total  strain  increment  to  total  stress  increment. 

The  important  difference  to  note  between  the  above  two  approaches  is  that  in  the 
initial  strain  method,  element  stiffness  matrices  are  computed  only  at  the  beginning  of  the 
analysis  whereas  for  the  tangent  stiffness  method,  they  have  to  be  computed  at  each  load 
increment.  Thus,  if  only  inelastic  material  is  being  studied  and  it  can  be  safely  assumed  that 
the  displacements  will  remain  small,  the  initial  strain  approach  is  found  to  be  much  more 
economical.  However,  if  nonlinear  geometric  behavior  must  be  incorporated  in  the  analysis, 
the  usual  procedure  has  been  to  recompute  element  stiffness  matrices  at  each  load  step.  Thus 
there  would  be  no  particular  advantage  in  using  the  initial  strain  approach.  However,  we 
shall  discuss  this  point  later. 

1 .2.2  Large  Deformations 

Mallet  and  Marcal  [4]  have  described  the  three  alternative  formulations  that  have  been 
used  in  the  past  for  solving  geometric  nonlinear  structural  problems  by  the  finite  element 
stiffness  method.  The  first  approach  seeks  a  displacement  state  which  minimizes  the  total 
potential  energy.  The  second  method  is  a  direct  one  where  the  equilibrium  equations  are 
solved  for  the  desired  displacements.  Since  these  equations  are  nonlinear,  such  techniques  as 
the  Newton -Raphson,  or  modified  Newton-Raphson,  method  have  been  used.  The  final  and 
most  popular  approach  is  the  linear  incremental  method.  The  equilibrium  equations  are 
linearized  in  each  increment  of  loading  to  obtain  stiffness  matrices  which  are  dependent  on 
the  current  state  of  the  structure.  This  is  the  approach  used  here  and  will  be  described  in 
detail. 

It  is  important  to  note  that  the  above  description  of  approaches  for  solving  nonlinear 
geometric  problems  is  the  first  line  classification  of  the  available  methods.  That  is,  the  three 
approaches  are  fundamental  and  any  solution  procedure  is  applied  to  the  total  potential 


*Numbers  in  brackets  designate  references  listed  on  page  94. 
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energy,  the  nonlinear  equilibrium  equations,  or  the  incremental  equations  obtained  by 
linearizing  the  equilibrium  equations. 

The  wide  popularity  of  the  linearized  incremental  formulation  is  apparently  due  to 
its  computational  advantage  over  the  other  two  formulations.  Efficient  linear  equation  solvers 
have  been  developed  and  can  be  easily  applied  to  the  linearized  formulation.  Moreover,  the 
incremental  formulation  allows  the  straightforward  inclusion  of  constitutive  relationships  which 
are  path  dependent. 

The  development  of  a  correct  linear  incremental  formulation  did  not  come  easily. 

Martin  surveyed  the  progress  up  to  1965  at  the  first  Wright-Patterson  meeting  on  matrix 
methods  [5]  and  provided  a  derivation  which  helped  to  clarify  matters.  Since  then,  the 
formulations  have  stabilized  and  the  correct  forms  are  understood  quite  well. 

The  early  derivations  of  the  incremental  formulation  were  based,  more  or  less,  on 
intuitive  approaches  where  equilibrium  was  examined  at  individual  nodal  points  by  using  local 
coordinate  frames  which  translated  and  rotated  with  the  deformation  of  the  structure.  The 
result  was  that  the  derivation  of  the  geometrical  or  initial  stress  matrix  for  the  same  element 
varied  from  author  to  author.  Martin  pointed  out  these  inconsistencies  and  derived  the 
correct  formulation  for  the  initial  stress  matrix  for  several  elements  by  using  an  initial  con¬ 
figuration  reference  frame.  However,  Martin  truncated  the  expressions,  dropping  terms  which 
proved  later  to  be  significant. 

If  the  Lagrangian  description  is  used  in  the  derivation,  as  was  done  by  Martin,  it  was 
determined  later  that  an  additional  matrix,  termed  by  Marcal  [6],  the  initial  displacement 
matrix  must  be  included  in  the  stiffness  equation.  Marcal  showed  the  importance  of  includ¬ 
ing  this  additional  matrix  and  presented  a  formulation  for  its  derivation. 

Hibbitt,  Rice,  and  Marcal  [7]  included  the  effects  of  following  loads  in  the  Lagrangian 
formulation  by  introducing  the  initial  load  matrix.  This  term  was  later  derived  by  other 
investigators,  e.g.,  Oden  and  Keys  [8],  but  apparently  its  effect  has  not  been  determined 
because  of  the  awkwardness  of  introducing  the  terms. 

The  most  recent  addition  to  the  incremental  stiffness  equation  has  been  the  equilibrium 
correction  term,  apparently  first  introduced  by  Hofmeister,  Greenbaum,  and  Evensen  [9].  A 
formulation  was  introduced  later  by  Stricklin  et  al.  [10].  It  will  be  shown  subsequently  that 
the  term  is  needed  because  the  pointwise  equilibrium  requirement  within  the  element  is  not 
satisfied.  The  addition  of  the  term  in  the  linearized  incremental  stiffness  equation  can  be 
viewed  as  including  the  zero1*1  term  in  a  Taylor  expansion  as  opposed  to  dropping  the  term  in 
a  simple  Euler  integration. 

Hibbitt  [11]  pointed  out  that  when  the  equilibrium  correction  term  is  included  in 
the  linear  formulation,  care  must  be  taken  to  ensure  that  the  incremental  stresses  are  summed 
properly.  He  showed  the  importance  of  using  the  expression  for  strain  (which  is  quadratic  in 
the  computed  incremental  displacements)  in  calculating  the  incremental  stresses.  Otherwise, 
the  equilibrium  correction  term  obtained  by  using  the  linearized  formulation  can  give  very 
poor  results. 
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1.3  GENERAL  SHELL  FINITE  ELEMENTS 


The  application  of  the  finite  element  method  to  shell  structures  occurred  early  after 
the  original  development  of  the  method;  see  Greene  et  al.  [12].  This  was  a  natural  applica¬ 
tion  since  solutions  for  a  large  number  of  practical  problems  could  be  obtained  only  by  some 
other  numerical  means  such  as  the  finite  difference  method.  These  early  applications  used 
Hat  plate  elements  to  model  the  curved  surface;  these  proved  to  be  adequate  if  one  was  will¬ 
ing  to  pay  the  price  for  inordinate  amounts  of  computing  time.  No  known  attempts  have  been 
made  to  apply  these  flat  elements  to  nonlinear  shell  problems,  and  that  is  just  as  well  since 
the  results  would  probably  have  been  disastrous. 

The  previous  application  of  flat  plate  elements  to  shell  problems  revealed  that  the  po¬ 
tential  for  solving  very  difficult  problems  was  there  [13,  14],  It  remained  to  improve  the 
method  so  that  it  could  be  applied  with  reduced  computing  time  and  cost  and  to  other  than 
linear  elastic  problems.  The  most  logical  improvement  needed  was  the  development  of  an 
adequate  curved  shell  element  with  reasonable  convergence  characteristics.  Unfortunately, 
this  proved  not  to  be  easy.  Although  elements  have  been  developed  that  are  certainly  more 
than  adequate,  there  is  still  no  one  curved  shell  element  that  is  accepted  universally. 

Gallagher  [15]  has  outlined  the  four  difficulties  that  have  impeded  the  derivation  of 
curved  shell  elements: 

a.  The  choice  of  an  appropriate  shell  theory  -  both  shallow  and  nonshallow  shell 
theories  have  been  used.  The  choice  of  a  particular  shell  theory  should  be  preceded  by  the 
decision  on  whether  to  compute  element  stiffnesses  in  local  element  coordinates  or  in  global 
coordinates.  A  direction  that  is  often  taken  is  to  use  shallow  shell  theory  and  then  compute 
individual  element  stiffnesses  in  global  coordinates,  that  is,  element  displacements  and  strains 
are  in  terms  of  global  coordinates.  This  limits  the  application,  then,  to  shells  that  are  reason¬ 
ably  shallow  and  thus  restricts  their  application.  On  the  other  hand,  nonshallow  shell  theory 
can  cause  other  problems  because  of  the  degrees  of  freedom  used  for  the  element.  These 
problems  arise  when  branched  shells  are  being  studied. 

b.  Description  of  the  geometry  of  the  element  —  depending  on  which  shell  theory  is 
used,  various  geometrical  quantities  are  required  in  the  strain -displacement  relationships. 
Curvatures  are  usually  required  but  not  always.  Few  investigators  have  recognized  the  im¬ 
portance  of  having  these  geometric  quantities  as  accurate  as  possible.  The  approaches  have 
been  varied,  but  one  feature  seems  common  to  most  developments  that  have  been  reported 
previously.  Along  with  assuming  a  certain  displacement  distribution  in  the  element,  the  sur¬ 
face  is  taken  to  be  some  polynomial  form  with  undetermined  coefficients.  The  values  of 
these  coefficients  are  determined  by  evaluating  geometric  quantities  within  the  element  and 
on  the  boundary.  The  form  of  the  polynomial  can  be  determined  uniquely  taking  a  sufficient 
number  of  these  points,  depending  on  the  order.  There  apparently  has  been  only  one  report 
of  an  investigation  being  conducted  to  determine  the  magnitude  of  errors  that  are  introduced 
by  taking  the  above  approach  [16]. 
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c.  Representation  of  rigid  body  modes  of  behavior  —  many  of  the  elements  that  have 
been  developed  to  date  are  seriously  deficient  in  that  the  assumed  displacement  patterns  do 
not  include  all  rigid  body  modes.  The  result  is  that  strains  are  induced  when  the  element 
undergoes  a  rigid  body  rotation  or  translation.  This  problem  is  usually  associated  with  ele¬ 
ments  that  are  represented  in  curvilinear  coordinates  where  rigid  body  motion  is  represented 
only  by  a  coupling  of  the  individual  displacements. 

d.  Satisfaction  of  requirements  related  to  continuity  of  displacements  —  to  achieve 
this,  the  u,  v,  and  w  displacements  must  be  represented  by  the  same  order  of  ploynomials. 

For  flexure  requirements,  the  first  derivative  of  the  displacements  across  element  boundaries 
must  be  continuous.  The  early  applications  of  flat  plate  elements  did  not  satisfy  either  of 
these  requirements  but  could  be  achieved  in  the  limit  as  the  grid  pattern  was  made  finer. 

There  has  been  an  extensive  development  of  curved  shell  elements  for  axisymmetric 
and  arbitrary  shells  [16-35].  The  various  approaches  can  best  be  described  by  discussing 
representative  examples;  the  remaining  elements  are  variations  on  these  basic  types. 

One  of  the  more  sophisticated  elements  developed  to  date  is  that  by  Argyris  and 
Scharph  [17].  They  used  fifth  order  polynomials  to  develop  a  triangular  element  with  1 8 
degrees  of  freedom  at  each  vertex  node  and  three  degrees  at  each  midside  node.  The  formu¬ 
lation  for  the  element  is  derived  in  triangular  surface  coordinates.  The  particular  form  of 
differential  geometry  is  derived  for  determining  the  metric  of  the  undeformed  and  deformed 
surface.  These  are  used  to  obtain  the  strain  displacement  relationships.  The  authors  report 
that  the  element  satisfies  continuity  requirements  and  that  no  strains  are  induced  in  the  ele¬ 
ment  under  rigid  body  motion. 

The  triangular  element  developed  by  Dupuis  [18],  based  on  a  previous  formulation  by 
Dupuis  and  Goel  [19],  has  been  used  for  a  wide  range  of  practical  applications.  The  element 
has  nine  degrees  of  freedom  at  each  node,  and  third  order  Lagrange  polynomials  are  used  to 
interpolate  the  displacements  within  the  element.  The  accurate  shell  theory  due  to  Koiter 
[36]  and  Budiansky  and  Sanders  [37]  is  used  for  the  strain  displacement  expressions.  This 
element  satisfies  the  continuity  and  rigid -body  strain  requirements. 

The  elements  that  probably  have  received  the  most  attention  are  those  developed  by 
Zienkiewicz  [20,  21].  The  elements  were  originally  reported  for  applications  to  thick  shells 
and  have  since  been  specialized  for  anlayses  of  thin  shells.  Both  triangular  and  quadrilateral 
shapes  have  been  developed.  The  elements  evolved  from  the  original  isoparametric  concept 
where  geometry  and  displacement  modes  are  assumed  to  have  the  same  functional  form.  The 
superparametric  concept  is  used  for  this  particular  application;  the  geometry  of  the  element  is 
defined  by  higher  order  polynomials  than  those  used  in  defining  the  displacements  within  the 
element.  No  shell  theory,  as  such,  is  used  in  developing  the  element.  The  basic  three- 
dimensional  state  of  stress  is  assumed  to  exist  throughout  the  element  except  that  the  normal 
stress  through  the  thickness  is  set  equal  to  zero  to  conform  to  shell  behavior.  There  are  five 
degrees  of  freedom  at  each  of  the  nodes  which  may  be  arbitrary  in  number,  depending  on  the 
order  of  the  polynomials  assumed  for  the  displacements.  The  rigid -body  condition  is 
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satisfied  by  the  element.  There  has  been  some  difficulty  in  applying  the  element  since  it  has 
a  tendency  to  be  too  stiff.  This  problem  has  been  circumvented  by  arbitrarily  reducing  the 
number  of  points  used  in  the  numerical  integration  procedure  for  computing  the  element  stiff¬ 
ness  matrix. 

Appropriate  shell  theory  has  been  employed  to  develop  elements  that  include  deforma¬ 
tion  due  to  transverse  shear.  This  approach  was  used  by  Key  and  Beisinger  [16]  to  develop  a 
quadrilateral  element  where  the  shell  theory  is  that  given  by  Washizu  [38],  The  advantage  of 
this  approach  is  that  the  continuity  requirements  demand  only  that  u,  v,  and  w  be  continuous. 
The  usual  shell  theory  requirement  that  normals  must  remain  normal  is  relaxed,  and  thus  the 
first  derivative  of  the  normal  displacement  w  does  not  have  to  be  continuous  across  element 
boundaries.  This  approach  makes  the  selection  of  displacement  interpolation  polynomials 
easier  and  Hermite  polynomials  are  used.  There  are  seven  degrees  of  freedom  at  each  of  the 
four  corner  nodes.  A  mapping  similar  to  that  due  to  Zienkiewicz  is  used,  and  the  element  is 
limited  to  shells  of  revolution.  It  does  not  satisfy  the  rigid -body  mode  requirement. 

The  approach  of  using  shallow  shell  theory  in  global  coordinates  is  given  by  Cowper 
et  al.  [22],  The  strain  displacement  expressions  of  Novozhilov  [39]  are  employed  in  the 
element  along  with  cubic  distribution  of  the  displacements  in  the  surface  of  the  shell  and 
fifth-order  polynomial  distribution  for  the  normal  displacement.  There  are  twelve  degrees 
of  freedom  at  each  node.  The  surface  of  the  shell  is  assumed  to  be  quadratic  in  form;  this 
dictates  that  the  element  must  have  constant  curvatures.  The  primary  disadvantage  of  this 
method  is  that  complete  shells  cannot  be  studied. 

A  variation  on  the  above  approach  was  also  presented  by  Cowper  et  al.  [23].  The 
element  was  assembled  in  a  global  system  on  the  shell  surface,  that  is,  the  element  stiffness 
was  first  computed  in  its  own  local  system  and  then  rotated  on  the  shell  surface  to  line  up 
with  the  Gaussian  coordinate  system.  This  is  different  from  the  first  developed  where  the 
elements  were  rotated  to  line  up  with  the  coordinate  system  in  the  projected  plane  of  the 
total  shell.  Thus  this  second  approach  removed  the  restriction  that  only  parts  of  complete 
shells  could  be  studied.  The  same  degrees  of  freedom  were  used  as  those  reported  for  the 
initial  Cowper  approach  [22], 

The  development  of  an  additional  triangular  shell  element  was  also  reported  by  Cowper 
et  al.  [23].  Few  details  were  given,  but  the  development  followed  previous  approaches  where 
the  degrees  of  freedom  and  stiffness  were  expressed  in  terms  of  the  curvilinear  coordinates  in 
the  shell  surface.  In  all,  36  degrees  of  freedom  are  used  in  the  element  with  12  at  each  of 
the  three  nodes.  The  authors  state  that  exact  rigid -body  modes  are  not  obtained  with  the 
element  and  that  this  point  has  often  been  overemphasized.  This  is  contrary  to  a  majority  of 
the  research  reported  previously;  see  Gallagher  [15]. 

The  method  of  using  shallow  shell  equations  in  local  element  coordinates  and  trans¬ 
forming  to  the  surface  was  also  reported  in  Strickland  and  Loden  [24]  where  a  triangular 
element  is  developed.  The  primary  difference  between  this  work  and  that  reported  by  Cowper 
et  al.  [23]  is  that  the  order  of  displacement  functions  is  less  and  thus  results  in  fewer  degrees 
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of  freedom  at  the  nodes.  Linear  functions  are  used  for  the  membrane  displacements  and  a  cubic 
function  for  the  normal  displacement.  There  are  a  total  of  15  degrees  of  freedom,  five  at  each 
node.  As  expected,  the  convergence  rate  to  the  correct  solution  for  increasing  number  of  ele¬ 
ments  is  less  than  that  reported  in  Cowper  et  al.  [22  and  23]. 


CHAPTER  2 

GENERAL  INCREMENTAL  FORMULATION 


We  turn  now  to  the  derivation  of  the  incremental  formulation  that  will  eventually  be 
implemented  by  the  finite  element  method.  Although  we  will  show  how  the  formulation  can 
be  applied  in  terms  of  finite  elements,  it  is  not  limited  to  this  method  but  could  be  imple¬ 
mented  by  another  numerical  approach  such  as  the  finite  difference  method.  The  derivation 
in  the  beginning  will  be  independent  both  of  the  magnitude  of  deflection  or  strain  and  of  the 
constitutive  relationship  between  stress  and  strain.  The  limitations  on  material  behavior  will 
be  described  later. 

2.1  DEFINITION  OF  STRAIN  AND  STRESS 

The  motion  of  a  point  in  a  body  is  described  by  a  one-to-one  mapping  of  the  type 

xi  =  X;  +  Uj  (2.1) 

where  Xj  is  the  initial  position  of  the  point  and  x;  is  its  position  at  a  later  time.  For  our  work 
here,  we  will  say  that  they  are  both  in  the  same  Cartesian  coordinate  system.  The  amount  that 
the  point  moves  is  u^  There  are  two  classical  methods  by  which  the  motion  is  described,  the 
Lagrangian,  or  material  description,  and  the  Eulerian,  or  spatial  description.  In  the  first, 
various  quantities  are  described  in  terms  of  X  as  the  independent  variable  whereas  in  the 
second,  x  is  the  independent  variable. 

For  the  measure  of  strain,  the  square  of  a  line  segment  before  and  after  deformation  is 
required.  These  are  given  respectively  by 

dS2  =  dXjdXj  =  5y  dXjdXj 

/dxi  \/9xj  \  (2.2) 

ds.,dx.dx.^_dx^_  dXkj 

where  5y  is  the  Kronecker  delta  and  the  summation  convention  is  used  with  the  range  from 
1  to  3  over  the  indices.  Thus  the  expression  for  strain  becomes 

ds2  -  dS2  =  dXjdXj  -  dXjdXj  =  2ey  dxjdxj  =  2Ey  dXjdXj  (2.3) 

The  symbol  e^  is  the  Eulerian  or  Almansi  strain  tensor  and  Ey  is  the  Lagrangian  or  Green  strain 
tensor. 

In  the  later  numerical  implementation  of  the  theory,  it  is  more  convenient  to  work  with 
the  material  description  and  thus  the  Green  strain  tensor  will  be  applied.  Using  equations 
(2.1),  (2.2),  and  (2.3),  the  expression  for  Green  strain  becomes 
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E,j  =  I/2(Uij  +  Uj_,  +  uk>1ukij) 


(2.4) 


where  the  comma  indicates  differentiation  with  respect  to  the  material  independent  variable  X.  i 

Since  the  derivation  and  application  of  the  methods  to  follow  will  be  performed  in  the 
material  coordinates,  it  will  be  necessary  to  have  a  stress  measure  which  is  in  terms  of  the 
undeformed  configuration.  This  will  be  the  Kirchoff,  or  Piola,  stress  tensor. 

We  consider  a  force  vector  df  acting  on  an  element  of  area  in  the  strained  configura¬ 
tion  and  want  to  assign  a  correspondence  rule  between  that  vector  and  one,  dF,  acting  on  the 
element  of  area  in  the  unstrained  configuration.  In  component  form,  we  write  the  Kirchoff 
rule  for  the  correspondence  [40]  as 


9Xi 

dFi=^dfJ 


(2.5) 


which  is  the  same  form  as  that  for  transforming  from  spatial  to  material  coordinates. 

Using  equation  (2.5),  we  can  obtain  the  result  that  we  actually  need,  i.e.,  the  relation¬ 
ship  between  true  or  Cauchy  stress  and  the  nominal  or  Kirchoff  stress.  The  derivation  is 
given  in  many  continuum  mechanics  text  books  and  will  not  be  repeated  here.  The  transfor¬ 
mation  of  Cauchy  stress  to  Kirchoff  stress  is 


8Xj  3X^ 
3xk  3xg  ° 


(2.6) 


where  is  the  Kirchoff  stress,  okg  is  the  Cauchy  stress,  p  is  the  density  in  the  deformed 
configuration,  and  p0  is  the  density  in  the  original  configuration.  The  Kirchoff  stress  tensor 
is  symmetric.  This  can  be  seen  above  since  the  Cauchy  stress  tensor  and  the  transformation 
are  symmetric. 


2.2  PRINCIPLE  OF  VIRTUAL  WORK  -  INCREMENTAL  THEORY 

The  method  under  consideration  for  solving  the  nonlinear  problems  will  have  to  be 
incremental,  i.e.,  obtained  in  finite  steps,  since  we  will  be  including  inelastic  material  effects 
which  are  path  dependent.  The  derivation  given  will  be  independent  of  material  properties 
but  can  be  applied  to  bodies  which  undergo  plastic  deformation,  as  will  be  described.  Al¬ 
though  we  will  eventually  be  considering  problems  that  are  characterized  by  small  strains,  the 
following  derivation  is  not  subject  to  that  limitation. 

The  three  configurations  in  the  deformation  path  that  are  of  interest  are  shown  in 
Figure  2.1:  the  original  undeformed  reference  configuration,  the  intermediate  deformed  con¬ 
figuration,  and  the  current  deformed  configuration  for  which  the  virtual  work  equation  will 
be  written.  The  intermediate  configuration  is  taken  as  the  one  just  before  the  current 
configuration. 
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CURRENT  DEFORMED 
CONFIGURATION 


Figure  2.1  —  Original  and  Deformed  Configurations 


Over  the  surface  of  the  body,  there  are  specified  displacements  and  surface  tractions,  F,  in 
terms  of  the  original  undeformed  reference  configuration.  There  are  also  body  forced,  P, 
expressed  in  terms  of  the  density  and  volume  of  the  undeformed  reference  configuration. 
The  virtual  work  equation  is  now  written  for  the  current  configuration 


f  (Sjj  +  ASij)5(AEij)dV  - J (p0Pi  +  p0APi)5(Aui)dV 
V 

-  (Fi  +  AFi  +  AGFi)6(Aui)dS  =  0 


(2.7) 


which  is  referenced  to  the  original  undeformed  state,  i.e.,  all  integrals  are  performed  over  that 
configuration  and  the  conditions  stated  above  with  respect  to  loads  hold  true.  The  symbols 
and  Ejj  have  been  defined  previously  and  are  the  Kirchoff  stress  tensor  and  the  Green  strain 
tensor,  respectively.  The  prefix  A  indicates  the  change  in  the  quantity  that  occurs  in  the  body 
in  going  from  the  intermediate  deformed  configuration  to  the  current  configuration.  The 
difference  in  the  surface  loading  caused  by  this  change  is  indicated  by  the  symbol  AqFj.  The 
subscript  G  means  that  the  change  in  force  is  associated  with  the  change  in  geometry. 

In  writing  equation  (2.7),  it  is  assumed  that  the  body  is  in  equilibrium  in  the  current 
configuration  under  the  body  forces,  the  external  applied  surface  tractions,  and  the  specified 
displacements.  The  body  is  then  assumed  to  execute  an  infinitesimal  virtual  displacement 
5(Au)  from  the  equilibrium  configuration  without  violating  the  specified  displacements.  It  is 
important  to  note  that  the  variation  is  only  on  the  displacements  with  the  prefix  A  since  the 
displacements  occurring  up  to  and  including  the  intermediate  deformed  configuration  just  be¬ 
fore  the  current  configuration  are  assumed  to  be  constants  in  the  equation. 

The  value  of  the  Green  strain  that  occurs  in  the  current  configuration  is 

Ejj  =  1/2  (Uj  +  AUj),j  +  1/2  (Uj  +  AUj),j  +  1/2  (um  +  Aum),j  (um  +  Aum),j  (2.8) 

with  the  variation  resulting  from  a  variation  in  AUj  given  by 

S(AEjj)  =  5 Ejj  =  1/2  5(AUj),j  +  1/2  5(AUj),j  +  1/2  um  j  5(Aum),j  +  1/2  umJ  6(Aum),j 

+  l/2(Aum),j  5(Aum),j  +  1/2  (Aum),j  6(Aum),j  (2'9) 

Now,  substituting  equation  (2.9)  into  equation  (2.7)  and  recognizing  the  symmetry  of 
the  strain  tensor,  we  may  write 


J (S^  +  ASjj) [  1  /2  (5(AUj),j  +  5( AUj),j)  +  um j  5(Aum  ),j  +  (Aum  ),j  5( Aum  )(j )  dV 
-/  (p0Pj  +  p0APj)5(Aui)dV-y'(Fj  +  AFj+  AGFj)5(Auj)dS  =  0 


(2.10) 


f 
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This  equation  can  now  be  rearranged  to  obtain 


r  (1)  (2) 

J  ASy[l/2  (5(Aui),j  +  5(Auj),i  +  um  i8(Aum),j]dV+y  Sij(AUm)>i6(Aum),j  dV 


+ 

v 


r  O)  r  (4)  ^  (5) 

j  ^ij(Aum),i5(Aum),jd  V-y  AGFi6(Aui)dS=J(p0Pi  +  p0APi)5(Aui)dV  (2  , , . 


r  (6)  .  (7) 

j  (Ft  + AFi)8(Aui)dS- j  S^l/2  ^(Au^  +  6(Auj),i)  +  um  i8(Aum),j]dV 


At  this  point  it  is  convenient  to  discuss  the  form  of  equation  (2.1 1).  We  will  substitute 
shortly  for  an  increment  of  stress  in  the  first  and  third  integrals  and  eliminate  the  virtual  dis¬ 
placements  when  the  above  equation  is  cast  into  the  form  for  the  finite  element  method.  The 
present  form  of  the  equation  is  such  that  we  may  compare  the  various  terms  with  the  linearized 
incremental  formulations  that  have  been  derived  by  other  investigators. 

Except  for  term  3,  each  of  the  above  numbered  terms  in  equation  (2.11)  is  an  equiva¬ 
lent  term  of  the  linearized  incremental  stiffness  equation  of  the  finite  element  method.  Term  1 
can  be  shown  to  be  the  small  displacement  plus  initial  displacement  matrix  due  to  Marcal  [6], 
Term  2  is  the  initial  stress  or  geometric  stiffness  matrix  [5].  Term  4  is  the  initial  load  matrix 
of  Hibbitt  et  al.  [7]  and  Oden  and  Keys  [8],  Terms  5  and  6  are  the  body  force  load  plus 
increment  of  load  and  the  surface  load  plus  increment  of  load,  respectively.  Term  7  is  the 
equilibrium  load  correction  term  of  Hofmeister  et  al.  [9]  and  Stricklin  et  al.  [10].  Apparently, 
term  3  has  never  been  included  and  applied  in  previous  incremental  nonlinear  finite  element 
analyses. 

When  we  substitute  for  an  increment  of  stress  in  equation  (2.1 1),  we  will  obtain  a  form 
of  the  incremental  equations  that  is  altogether  different  from  the  linear  form  that  has  been 
used  by  previous  investigators.  It  is  primarily  term  3  although  that  is  lost  in  the  linearization 
process. 

Before  we  go  on,  note  one  further  point  concerning  equation  (2.1 1).  The  equilibrium 
load  correction  term,  term  7,  can  be  written  as 

•  +  ^ij  um,i  ^  9 

v  v 

but 


xm,i  =  5im  +  um,i 
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So  term  7  becomes 


/S ij  xm,i5(Aum)>jdV 

v 

Now  consider 

/ (S,j  xm,i  6(Aum)),j  dV  =/( Sy  8(Aum)dV+y Sjj  xm 8(Aum),j  dV 

V  V 

and  by  the  Gauss  integral  theorem 

/lS ij  xm,i  5(Aum)l  -j  dV=/ Sij  Vi  ^6(Aum)dS 

V  s 

where  nj  is  the  component  of  the  unit  normal  to  the  surface.  Thus,  term  7  may  be  written 
/Sij  xm,i^Aum)’j  xmi),j  5(Aum)+^Sjj  xrn  ir,j5(Aum)dS 

V  vs 

We  now  take  term  7  and  add  it  to  the  first  part  of  term  5  and  the  first  part  of  term  6; 
regrouping  and  paying  particular  attention  to  the  signs,  we  can  write  the  expression 

/ KSy  V i>'j  +  Po Pm  1  8< Aum) d v  +/( - Sy  xm,i>'j  +  Fm)6(ium)dS 
V 

At  each  step  of  the  loading,  we  assume  a  displacement  pattern  which  results  in  an  assumed 
strain  and  stress  distribution.  If  this  resulting  stress  pattern  satisfied  the  equations  of  equilib¬ 
rium  at  each  point  in  the  body  and  at  the  surface,  we  could  write  the  following  equations 
(expressed  in  the  Lagrangian  description) 

^ij  xm,  i^’j  +  ~  ® 

(2.12) 

^ij  xm ,  i  ~^m 

The  point  is  that  these  conditions  are  not,  in  general,  satisfied  and  that  is  why  the  above- 
considered  terms  have  to  be  included  in  equation  (2.11).  This  close  examination  of  the  virtual 
work  equation  points  out  the  necessity  of  including  all  of  the  terms  considered.  It  is  often 
stated  that  equilibrium  is  not  satisfied  at  any  given  load  level  in  the  incremental  solution  to 
nonlinear  problems.  Specifically,  it  is  the  pointwise  equilibrium  requirements  of  equation 
(2.12)  that  are  not  being  satisfied. 

We  now  propose  to  rearrange  equation  (2.11)  into  a  form  which  will  allow  a  different 
solution  method  to  nonlinear  incremental  problems.  This  is  shown  as  follows 
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f  ASjj  [1/2  (SCAuj)^  +  5( AUj ),j)  +  (Um  i  +  Aum  j)  5(  Aum  ),j  ]  d V 
v 

Pi  +  p0APj)8(Aui)dV-»y(Fi  +  AFi+AGFi)5(Aui)dS  (2  ,3) 

s  s 

-f  Sy  [1/2  (6(AUj),j  +  6(Auj),j)  +  (um  i  +  Aum  j)8(Aum),j]  dV 
v 

The  first  difference  to  be  noted  between  equations  (2.1 1)  and  (2.13)  is  that  term  3  of  equation 
(2.1 1)  is  put  into  the  first  term  of  that  equation.  Secondly,  the  initial  stress  term,  term  2,  of 
equation  (2.1 1)  is  moved  to  the  right  side  and  put  into  term  7.  Finally,  the  initial  load  term 
of  equation  (2.1 1)  is  put  into  term  6  of  that  equation. 

The  initial  load  term  of  equation  (2.1 1)  is  a  function  of  the  change  in  area  which  is  a 
function  of  the  change  in  displacements  that  occur  in  each  increment.  When  it  is  kept  on  the 
left  side  and  its  contribution  added  to  the  stiffness  matrix,  the  matrix  is  nonsymmetric  and 
causes  problems  in  the  solution  process.  Thus  it  is  more  convenient  to  move  the  term  to  the 
right  side  and  include  its  contribution  as  part  of  the  iteration  on  displacements  that  is  per¬ 
formed  to  reach  convergence.  This  is  suggested  in  Hibbitt  [11], 

We  now  complete  the  formulation  of  equation  (2.13)  by  first  expressing  an  increment 
of  stress  in  terms  of  an  increment  of  strain  by  the  constitutive  relationship 

ASij  =  DijabAEab  (2.14) 

where  Dyab  is  a  known  function  of  the  current  state.  As  described  previously,  both  the  stress 
and  strain  are  measured  with  respect  to  the  undeformed  configuration.  In  writing  equation 
(2.14),  no  limitation  is  put  on  the  magnitude  of  strain.  We  ask  only  that  the  relationship  be¬ 
tween  stress  and  strain  remain  linear.  If  the  material  goes  plastic,  then  we  will  require  the 
strain  to  be  small  and  will  use  the  constitutive  relationship  derived  in  the  next  chapter. 

The  increment  of  strain  can  be  computed  by 

AEab  =  1 12  (ua  +  Aua ),b  +  1/2  (ub  +  Aub  ),a  +  1  /2  (um  +  Aum  ),a  (um  +  Aum  ),b 
-|/2(ua)b+ub)a)-l/2umauni>b 


where  the  increment  of  strain  is  that  in  the  current  configuration  minus  the  strain  in  the 
intermediate  state  just  before  the  current  state.  The  above  reduces  to 


AEab  =  1/2  (Aua  b  +  Aub  a)  +  1/2  um  a  Aum>b  +  1  /2  Aum>a  um  b  +  1/2  Au 


m,a  ^um,b 


(2.15) 


After  recognizing  the  symmetry  in  a  and  b,  the  expression  for  an  increment  of  stress 
now  becomes 


ASij  =  Dijab  l 1  / 2  (Aua,b  +  Aub,a>  +  um,a  Aum,b  +  X>2  Aum,a  Aum,b  1  <2- 1 6) 

Now  substituting  equation  (2.16)  into  the  left  side  of  equation  (2.13),  we  obtain 

/•  (1) 

J  Dijabf  !/2  <Aua,b  +  Aub,a)  +  ( Va  +  !/2  Aum,a)Aum,b  1 

V 

[  1  /2  (5(Aui),j  +  6(Auj),i)  +  (un>i  +  Aun  i)  6(Aun),j  ] 

f  f  ^  (2  17) 

= J  (p0  Pj  +  p0  APj)5(AUj)dV  +J  (Fj  +  AFj  +  AG Fj) 5( AUj)dS  U’  1 

V  s 

r  (4) 

-ySij[l/2(6(Aui),j+6(Auj),i)  +  (umii  +  Aumi)6(Aum),j]dV 


which  is  the  final  form  of  the  virtual  work  equation.  Note  that  the  left  side  of  the  equation, 
which  corresponds  to  the  stiffness  matrix,  is  nonsymmetric.  Since  we  will  want  to  use  the 
efficiency  of  a  symmetric  equation  solver,  it  will  be  convenient  to  separate  the  term  into  a 
symmetric  part  and  a  nonsymmetric  part.  This  will  be  described  subsequently. 


2.3  FINITE  ELEMENT  FORMULATION  OF  THE  EQUILIBRIUM  EQUATION 

Equation  (2.7)  was  written  without  mentioning  the  origin  and  orientation  of  the 
material  reference  frame.  It  is  often  convenient  to  have  a  reference  frame  for  each  individual 
element  when  the  finite  element  method  is  used.  Thus,  the  previous  derivation  could  be 
thought  of  as  being  done  for  a  single  element  with  its  own  coordinate  system.  Eventually  we 
will  want  to  be  able  to  assemble  all  of  the  elements  into  a  single  system.  This  causes  no 
problem  since  all  of  the  quantities  in  equation  (2.7)  can  be  transformed  easily  because  they 
are  tensor  quantities.  Thus  we  will  not  make  a  distinction  between  the  local  element  system 
and  the  common  global  system.  This  is  necessary  only  when  we  start  to  apply  the  method  to 
a  particular  problem. 

The  finite  element  stiffness  approach  is  based  on  the  approximation  of  representing  the 
displacment  at  any  point  in  an  element  in  terms  of  polynomials.  For  our  application  here,  the 
polynomials  will  be  written  in  terms  of  the  material  coordinates.  The  coefficients  in  the  poly¬ 
nomials  are  represented  as  linear  combinations  of  the  displacements  that  occur  at  the  nodal 
points.  Thus  we  seek  a  relationship  where  the  displacement  at  any  point  in  the  element  can 
be  determined  in  terms  of  the  displacements  that  occur  at  the  nodes. 

To  determine  the  relationship  between  the  displacements  at  the  nodes  and  interior 
points,  we  start  by  selecting  the  polynomial  form  that  will  be  used  for  the  three  u,  v,  and  w 
displacements  at  the  interior  points.  We  represent  this  by 

ui  =  rijaj  (2.18) 
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where  is  an  array  of  polynomial  terms  and  aj  are  unknown  coefficients.  The  range  on  i  is 
one  to  three  and  the  range  on  j  is  from  one  to  the  number  of  nodal  points  in  the  element 
times  the  number  of  displacements,  i.e.,  degrees  of  freedom  at  each  nodal  point. 

Next  a  set  of  generalized  displacements  is  introduced  [41]  such  that: 

a.  The  connection  of  one  element  to  another  element  along  a  common  boundary  is 
expressed  by  equating  each  of  the  respective  generalized  displacements  pertaining  to  the 
boundary. 

b.  The  displacement  field  along  the  common  boundary  is  uniquely  defined  in  terms  of 
the  generalized  displacements  that  occur  at  the  nodal  points  on  the  boundary  and  thus  ensures 
the  continuity  of  displacements  across  the  boundary. 

c.  The  number  of  generalized  displacements  is  equal  to  the  number  of  unknown 
coefficients  in  equation  (2.18). 

The  set  of  generalized  displacements  is  expressed  by 


The  set  of  generalized  displacements  Vj  is  composed  of  the  displacements  u;  and,  in  general, 
their  derivatives.  Thus,  the  terms  of  the  array  Rjj  are  composed  of  the  terms  of  r^  plus  terms 
created  by  taking  the  derivatives  of  the  terms  with  respect  to  the  material  coordinates.  The 
introduction  of  the  generalized  displacements  sets  the  form  of  generalized  forces  which  evolve 
in  the  process  when  the  consistent  loading  vectors  are  derived. 

An  expression  is  required  for  the  coefficients  aj  in  terms  of  the  generalized  displace¬ 
ments  at  the  nodes.  This  is  obtained  by  first  substituting  the  coordinates  of  the  nodal  points 
into  equation  (2.19)  to  obtain  a  set  of  equations  equal  in  number  to  the  number  of  coeffi¬ 
cients  aj : 

Vj  =  Rjj  ^  (2.20) 

where  the  bar  indicates  that  the  equations  are  written  at  the  nodal  points.  This  set  of  equa¬ 
tions  is  then  solved  to  obtain 


a-  =  R'1  v-  =  A  -  V' 

i  ij  j  ij  j 


(2.21) 


and  is  the  desired  result  which  expresses  the  coefficients  of  the  polynomials  of  equation  (2.18) 
as  linear  combinations  of  the  generalized  displacements  that  occur  at  the  nodes.  Each  term  in 
Ay  is  a  constant.  Substituting  (2.21)  into  (2.18),  we  obtain 


ui  =  rijAjkvk 


(2.22) 


the  displacements  u,  v,  and  w  at  any  point  in  the  element  expressed  in  terms  of  the  generalized 
displacements  at  the  nodes. 
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Using  equation  (2.22),  we  may  now  obtain  the  finite  element  stiffness  equation  from 
the  virtual  work  equation  (2.17).  But  first  we  need  the  following  expressions 

ui,£  rij,£Ajkvk 

(2.23) 

Aui,£  ~  rij,£  Ajk  Avk 

(2.24) 

5(Aua)  =  rij,£Ajk5(A7k) 

(2.25) 

each  of  which  is  obtained  from  (2.22).  Equation  (2.25)  is  the  result  of  r^  being  a  function  of 
the  initial  position  only. 

Substituting  (2.23)  through  (2.25)  into  (2.17),  the  following  is  obtained  for  each  of 
the  terms  of  that  equation 

Term  1 

J*® ijabt  (rae,b  +rbe,a-^  +  ^d  +  ^  A''d^rhc,arhe,b^cdl 

V 

H/2  (rin  j  +  rjn  i)  +  (ve  +  Avg)rmk  irmn  jAk£]  dV  Anq  AefA?f5(Avq) 

(2.26) 

Term  2 

/ (PoPi  +  Co^i^ijiVAjqSCA^) 

V 

(2.27) 

Term  3 

/ (F,  +  AF,  +  AG  Fj)  Tjj  dS  Ajq  8(  Av^) 
s 

(2.28) 

Term  4 

Aj  1 1  /2  <rin,j  +  +  (*£  +  A^rmk,irmn,j  Akfil  dV  Anq  5(  A  V 

(2.29) 

V 


The  virtual  displacements,  8  (Avq),  of  the  nodal  points  are  arbitrary  and  thus  can  be  eliminated. 
In  other  words,  since  they  are  arbitrary,  an  equation  in  just  the  coefficients  of  the  virtual  dis¬ 
placements  must  hold  true.  The  result  of  equating  the  coefficients  is  the  stiffness  equation. 

An  important  point  to  remember  at  this  stage  is  that  the  virtual  displacement  of  a 
nodal  point  generalized  displacement  may  be  thought  of  as  causing  virtual  displacements  with¬ 
in  all  of  the  elements  that  are  connected  to  the  nodal  point.  Thus  all  of  these  elements  are 
represented  in  the  stiffness  equation  that  is  written  to  eliminate  the  generalized  virtual 


displacements  at  the  nodes.  We  noted  earlier  that  the  equations  may  be  transformed  easily  - 
so  there  is  no  problem  in  writing  an  equation  in  a  single  coordinate  system  in  the  coefficients 
of  the  virtual  displacements. 

The  matrix  form  of  the  stiffness  equation  can  finally  be  obtained  by  letting  the  first 
bracketed  part  in  (2.26)  be  represented  by  the  matrix  [Bl]  and  the  second  bracketed  part  by 
the  matrix  [B2].  The  resulting  matrix  form  is 


([Aiy  [B2]T[D][Bl]dV[A])j  Av}=[A] 

v  v 

+  [A]y'[r]T(  |F|  +  |AF}  +  |AGF[)dS-[A]y'[B2]T  |sf  dV 


y*[r]T(p0jp}  +  p0  |AP})dV 


s 


V 


(2.30) 


where  the  stress  in  (2.29)  has  been  written  as  a  column  vector.  (The  change  in  surface  trac¬ 
tions  due  to  a  change  in  geometry  will  be  derived  later  for  the  quadrilateral  element.)  The 
first  two  terms  on  the  right  are  the  generalized  forces;  their  form  is  set  by  the  choice  of 
generalized  displacements. 

In  forming  the  stiffness  matrix  and  the  equilibrium  correction  term  (the  first  and  last 
terms  of  (2.30),  respectively),  it  is  noted  from  (2.26)  and  (2.29)  that  the  increment  of  dis¬ 
placement  that  occurs  in  going  from  the  intermediate  deformed  configuration  to  the  current 
configuration  is  required.  Again,  these  are  the  quantities  denoted  by  Av  .  Since  the  increment 
of  displacement  is  the  fundamental  unknown,  some  method  is  required  to  estimate  the  values 
to  be  used  in  forming  each  of  the  terms.  This  estimate  is  obtained  by  a  linear  extrapolation  of 
the  increment  of  displacements  that  occur  in  the  previous  increment.  Thus,  if  the  load  steps 
are  equal,  for  example,  then  the  estimate  of  the  displacements  is  set  equal  to  those  that 
occurred  in  the  previous  increment.  An  iteration  is  carried  out  at  each  load  step  to  reach  con¬ 
vergence  (within  a  given  tolerance);  the  details  will  be  described  later. 

As  noted  previously,  the  stiffness  matrix  derived  above  is  nonsymmetric.  To  circum¬ 
vent  the  problem  of  having  to  solve  a  nonsymmetric  set  of  equations  after  the  elements  are 
assembled  into  the  common  global  system,  the  stiffness  matrix  is  decomposed  into  a  symmetric 
and  skew  symmetric  part,  that  is, 

Kij  =  l/2(Kij  +  Kji)+l/2(Kij-Kji)  (2.31) 

where  the  first  term  on  the  right  is  the  symmetric  part  and  the  second  is  the  antisymmetric 
part.  Thus,  simplifying  the  notation  of  equation  (2.30)  and  putting  the  nonsymmetric  part  of 
the  stiffness  matrix  on  the  right  side,  we  may  write 

[K] s  {Av}  =  |p}  +  { AP}  +  {f}  +  Iaf}  +  |  AGF}  -  |e[-  [K]  A  { Avest}  (2.32) 


The  subscripts  S  and  A  respectively  indicate  the  symmetric  and  antisymmetric  parts  of  the 
stiffness  matrix.  The  symbol  E  stands  for  the  equilibrium  correction  term  and  the  remaining 
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terms  are  self-explanatory.  The  nonsymmetric  part  of  the  stiffness  matrix  is  now  taken  as  an 
equivalent  set  of  forces  by  using  the  estimated  displacements  that  occur  in  going  from  the 
intermediate  configuration  to  the  current  configuration. 

2.4  INITIAL  STRESS  PROBLEMS 

We  close  this  chapter  with  a  brief  account  of  how  the  above  derivation  stands  with 
respect  to  previous  derivations;  the  closest  are  those  that  evolve  from  considering  the  incre¬ 
mental  problem  as  an  initial  stress  problem. 

The  initial  stress  problem  has  a  long  history  that  dates  back  to  Cauchy  [42].  More 
recently,  Green,  Rivlin,  and  Shield  [43]  and  Murnaghan  [44]  have  studied  the  problem.  The 
formulation  is  presented  in  the  textbook  by  Washizu  [38]. 

Finite  element  investigators  who  have  used  the  initial  stress  approach  include  Felippa 
[45],  Hofmeister,  Greenbaum,  and  Evensen  [9],  Murray  and  Wilson  [46],  Bergan  [47],  and 
Yaghmai  [48]. 

Briefly,  the  approach  that  seems  to  be  common  to  all  of  the  above,  except  that  of 
Hofmeister  et  al.  [9],  is  to  write  the  virtual  work  equation  for  two  successive  configurations, 
take  the  difference  in  the  two,  and  then  linearize  the  result  to  obtain  incremental  stiffness 
equations.  Various  coordinate  systems  are  used  in  the  derivations. 

The  approach  of  Hofmeister  et  al.  [9]  is  to  use  the  intermediate  deformed  configura¬ 
tion  shown  in  Figure  2.1  as  the  current  reference  frame  and  to  consider  this  as  the  initial  state 
with  the  body  subject  to  initial  stress.  Thus,  the  method  involves  updating  coordinates  and 
also  transforming  stresses  to  a  common  system  so  that  they  may  be  added.  The  importance 
of  this  last  point  was  brought  out  by  Nemat -Nasser  and  Shatoff  [49]  who  derived  an  incre¬ 
mental  approach  and  implemented  it  with  the  finite  difference  method. 

The  two  primary  points  to  be  said  for  the  derivation  presented  here  are,  first,  its 
simplicity  and,  second,  the  nontruncation  of  terms.  A  linearization  process  is  imposed,  as 
discussed  in  the  previous  section,  but  the  equations  are  complete  and  answers  can  be  obtained 
as  accurately  as  desired  by  using  iteration.  Moreover,  as  will  be  shown  later,  even  without 
iteration,  the  formulation  seems  to  give  better  results  than  have  been  obtained  by  previous 
incremental  formulations. 


CHAPTER  3 

CONSTITUTIVE  LAWS  OF  PLASTICITY 


This  chapter  reviews  the  basic  theory  of  plasticity  required  for  the  stress-strain  rela¬ 
tionship.  Many  publications  may  be  consulted  for  a  more  extensive  treatment,  e.g.,  the  text 
by  Hill  [50],  the  work  of  Koiter  [51],  and  the  extensive  review  by  Naghdi  [52].  Here  the 
description  will  consider  only  small  strain  problems  which  may  be  assumed  to  be  independent 
of  temperature  and  time. 

3.1  BASIC  ASSUMPTIONS  -  INITIAL  YIELD  CONDITION 

The  state  of  stress  is  discussed  in  terms  of  a  general  nine-dimensional  stress  space.  In 
this  space  there  is  an  envelope  which  defines  the  limit  of  elastic  behavior  of  a  point  within  a 
body.  As  load  is  slowly  applied,  the  state  of  stress  at  a  single  point  in  the  body  is  defined  by 
a  trajectory  which  moves  in  the  stress  space.  There  is  such  a  trajectory  for  all  points  in  the 
body.  As  long  as  the  trace  of  the  trajectory  moves  within  the  envelope  and  does  not  touch 
the  surface,  the  material  at  the  point  under  consideration  does  not  undergo  any  plastic  defor¬ 
mation.  Only  trajectories  that  are  within  and  on  the  boundary  of  the  envelope  have  any 
meaning. 

The  envelope  is  referred  to  as  the  initial  yield  surface;  it  can  be  shown  to  be  convex 
and  contains  the  origin  of  the  stress  space.  On  the  basis  of  experimental  data  and  the  be¬ 
havior  of  inelastic  deformations,  the  shape  of  the  surface  is  considered  to  be  limited  to  several 
general  configurations.  The  initial  state  of  the  material  is  usually  assumed  to  be  isotropic  and 
shows  no  preferred  direction  with  respect  to  stress.  Further,  it  is  known  that  the  initial  yield 
is  primarily  brought  about  by  shearing  stresses.  Bridgman  [53],  among  others,  has  shown 
experimentally  that  the  yielding  is  primarily  independent  of  hydrostatic  pressure.  It  is  there¬ 
fore  convenient  to  express  the  envelope  or  yield  function  in  terms  of  the  invariants  of  the 
deviatoric  stress  tensor 

aij  =  aij  ~  1/3  ^ij  wkk  (3.1) 

where  the  invariants  are 

Ji  =crii  =  0 

h  =  (3.2) 

J3  =  J/3  °!j  °jk  aki 

Thus,  the  initial  yield  function  may  then  be  expressed  as 


f=f(J2,J3) 


(3.3) 
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Many  functions  have  been  proposed  for  the  initial  yield  condition  but  two  have  become 
the  most  reasonable  for  simulating  experimental  data.  The  first  is  the  Tresca  function 
expressed  in  terms  of  maximum  shearing  stresses.  It  is  assumed  that  plastic  flow  can  occur  at 
the  point  where  the  maximum  shearing  stress  reaches  the  yield  stress  obtained  in  pure  shear. 
The  second  yield  function  is  that  due  to  von  Mises  and  is  the  one  that  has  been  applied  in  the 
work  here.  In  formulating  the  function,  von  Mises  ignored  the  effects  of  the  invariant  J3  and 
expressed  the  condition  by  the  equation 

f  ( J2  )  =  J2  -  1/3  ag  =  1/2  a[j  ojj  -  1  /3  a\  =  0  (3.4) 

where  oQ  is  the  yield  stress  obtained  in  a  simple  tension  test.  A  convenient  geometrical  inter¬ 
pretation  of  the  function  is  obtained  when  it  is  plotted  in  principal  stress  space.  It  turns  out 
that  the  yield  surface,  or  envelope,  is  a  circular  cylinder  with  the  axis  equally  inclined  to  the 
three  principal  stress  axes.  The  cylinder  has  infinite  length. 

3.2  LOADING  CRITERIA 

Once  the  state  of  stress  at  a  point  is  such  that  plastic  flow  can  occur,  a  criterion  must 
be  established  to  determine  whether  it  actually  takes  place.  Assume  now  that  the  yield 
function  for  the  general  case  becomes  a  function  of  the  plastic  strain  e^,  that  is, 

f=f(aij)efj)  =  °  (3.5) 

and  the  function  no  longer  defines  just  initial  yield  but  also  includes  additional  information 
and  is  now  called  the  loading  function. 

Consider  the  differential 


df  =  r~  do  -  +  dej*:  (3.6) 

3aij  1J  9efj  1J 

at  a  point  where  equation  (3.5)  is  satisfied.  Then  df<0  implies  that  f<0  and  is  the  condition 
for  unloading  from  a  plastic  state  to  an  elastic  state.  An  additional  requirement  for  unloading 
is  that  no  change  can  occur  in  the  plastic  strain.  The  terms  neutral  loading  and  loading  can  be 
defined  by  using  a  similar  argument.  In  summary, 


During  unloading: 

^djij<0’f=0 

(3.7) 

During  neutral  loading: 

df 

-d„ir°,f=o 

(3.8) 

During  loading: 

0 

11 

o' 

A 

•0 

(3.9) 
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A  geometric  interpretation  can  be  associated  with  each  of  the  above  since  the  loading  surface 
is  closed.  For  unloading,  the  stress  increment  vector  is  directed  inward;  during  neutral  loading, 
it  is  tangent  to  the  surface,  and  during  loading,  it  is  pointed  outward. 

3.3  WORK  HARDENING  -  THE  DRUCKER  POSTULATE 

When  a  uniaxial  test  specimen  exhibits  work  hardening,  the  stress  is  a  monotonically 
increasing  function  of  increasing  strain.  A  generalization  to  combined  stress  states  is  required. 
This  generalization  is  attributed  to  Drucker  and  is  fundamental  in  the  theory  of  plasticity 
since  all  of  the  information  required  for  the  stress-strain  relationships  can  be  derived  by  using 
the  concept.  Drucker  stated  that  if  a  body  in  equilibrium  has  its  stress  state  changed  con¬ 
tinuously  by  an  external  agency  such  that  the  final  stress  state  coincides  with  the  original 
one  (i.e.,  a  closed  cycle),  then  the  work  done  by  the  external  agency  is  not  negative. 

Assume  that  a  body  is  in  an  equilibrium  state  of  stress  a y  and  strain  ey.  Then,  if 
small  surface  tractions  are  applied  to  the  body,  small  increments  in  stress,  day,  and  strain, 
dey,  will  occur.  Once  the  loading  is  removed,  all  of  the  elastic  strains,  dey,  will  be  recovered. 
Then  the  body  has  experienced  work  hardening  if  the  following  two  conditions  are  true: 
during  loading; 

daijdeij>0  (3.10) 

and  on  completion  of  the  cycle, 

daij(deij-de?j)>0  (3.11) 

The  quantity  in  parentheses  is  the  nonrecoverable  plastic  strain  dey  and  thus  the  second  equa¬ 
tion  can  be  rewritten 

doydefpO  (3.12) 

3.4  FLOW  RULE  FOR  WORK-HARDENING  MATERIALS 

As  a  consequence  of  the  Drucker  postulate,  the  following  three  properties  of  the  load¬ 
ing  function  and  inelastic  deformations  can  be  obtained. 

a.  The  initial  yield  surface  and  all  following  loading  surfaces  are  convex. 

b.  At  a  regular  (smooth)  point  on  the  loading  surface,  the  strain  increment  vector 
must  be  directed  normal  to  the  surface. 

c.  An  increment  of  plastic  strain  must  be  a  linear  function  of  an  increment  of  stress. 
Statement  a  can  be  shown  to  be  true  by  considering  equation  (3.12).  If  each  of  the 

increments  of  stress  and  strain  is  associated  with  a  respective  vector,  then  equation  (3.12) 
states  that  the  dot  product  of  the  two  vectors  must  be  equal  to,  or  greater  than,  zero.  This 
means  that  the  angle  between  the  two  vectors  is  less  than,  or  equal  to,  90  degrees.  In  careful 
consideration  of  this,  it  may  be  simply  stated  that  there  can  be  no  inward  dents  on  the  load¬ 
ing  surface  and  thus  must  be  convex  at  all  points. 
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Statement  b  may  be  verified  by  using  an  argument  similar  to  that  given  above.  The 
additional  conclusion  can  be  reached  that  the  direction  of  the  strain  increment  vector  is 
independent  of  the  direction  of  the  stress  increment  vector.  Statement  b  may  be  expressed 
mathematically  by 

de^  =  A  r—  (3.13) 

1J  9ffij 

where  A  is  a  positive  function  dependent  on  stress,  strain,  and  strain  history.  It  is  noted  that 
the  function  f  may  be  viewed  as  a  plastic  potential. 

The  meaning  of  statement  c  can  be  shown  by  considering  the  following.  The  loading 
function  must  remain  zero  during  loading,  i.e.,  in  passing  from  one  plastic  state  to  another 
plastic  state.  This  is  expressed  by 

ffajj.  ejj)  =  0  (3.14) 

and  is  what  Prager  [54]  called  the  consistency  condition.  Thus,  we  must  have 

df=^-doii  +  —  de£  =  0  (3.15) 

9tJij  1J  3ejj  1J 

Now,  substituting  (3.13)  into  (3.15)  and  solving  for  A,  we  obtain 

9f  ^ 

9akC  °k ® 

A *  ~  df  a7~  (3  I6) 

aefj  a°ij 

Next  the  linear  relationship  between  an  increment  of  stress  and  an  increment  of  strain  is  ob¬ 
tained  by  substituting  (3.16)  into  (3.13): 


de 


p 

mn 


3f  3f 

9okC9ffmn  , 

3f  3f  d(7kC 

aef  a«ij 


(3.17) 


The  relationship  is  linear  since  the  fractional  quantity  does  not  contain  the  stress  increment 
vector.  This  characteristic  will  be  shown  subsequently  in  an  evaluation  of  the  factor. 

All  of  the  above  could  be  derived  by  the  more  general  theory  developed  by  Koiter 
[51  ]  for  loading  surfaces  with  singular  points,  i.e.,  pointed  surfaces.  This  generalization  is  not 
required  here  since  the  von  Mises  expression  is  used  for  the  initial  yield  surface  and  subsequent 
loading  surfaces. 
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3.5  SUBSEQUENT  LOADING  SURFACES  -  HARDENING  RULES 


It  is  necessary  to  know  how  the  loading  surface  changes  as  load  is  continuously  applied 
to  the  structure.  Two  principal  theories  have  been  used  to  determine  this  change,  the  isotropic 
rule  and  the  kinematic  hardening  rule  first  suggested  by  Prager  [55]. 

The  basic  concept  behind  the  isotropic  rule  is  that  there  is  a  continuous  growth  of  the 
loading  surface  as  plastic  loading  occurs  at  a  point  in  the  body  under  study.  This  is  in  sharp 
contrast  to  experimental  evidence  since  the  theory  predicts  a  negative  Baushinger  effect.  That 
is,  as  a  point  in  a  uniaxial  test  specimen  is  loaded,  say  into  the  inelastic  compression  range, 
and  then  unloaded  and  then  carried  into  the  tension  range,  the  point  at  which  plastic  flow 
starts  to  occur  is  at  a  greater  tensile  stress  than  if  the  specimen  were  first  loaded  into  the 
tension  range.  The  equivalent  concept  is  shown  in  Figure  3.1  in  principal  stress  space.  The 
important  point  to  note  is  that  the  origin  of  the  loading  surface  remains  fixed  and  the  surface 
grows  with  increasing  loading. 

The  kinematic  rule  is  more  reasonable  since  it  is  in  better  agreement  with  experimental 
data.  The  principal  idea  behind  this  concept  is  shown  in  Figure  3.2.  Both  stresses  and  strains 
are  associated  with  a  displacement  (reason  for  term -kinematic)  such  that  the  origin  of  the 
loading  surface  is  allowed  to  translate  while  the  size  and  shape  of  the  surface  is  not  allowed 
to  change.  The  motion  of  the  surface  must  be  considered  in  a  general  stress  space  with  con¬ 
straints  being  applied  for  special  conditions,  for  an  example,  no  motion  would  be  allowed  in 
the  third  zero  principal  stress  direction  for  plane  stress. 

3.6  EQUIVALENT  STRESS  AND  EQUIVALENT  STRAIN 

The  hardening  rules  for  combined  stress  are  related  to  the  simple  uniaxial  test  specimen. 
Some  measure  of  stress  and  strain  must  be  used  to  make  the  relationship  meaningful.  In  other 
words,  given  a  state  of  complex  stress  and  strain  at  a  point,  what  is  the  relationship  between 
this  state  and  our  experimental  data  for  a  uniaxial  test  specimen? 

A  definition  for  equivalent  stress  is  obtained  from  equation  (3.4).  This  is  the  require¬ 
ment  for  first  yielding  to  occur  under  a  condition  of  combined  stress.  Thus  the  effective 
stress  is 

ae  =  [3/2  a[j  o[-]1/2  (3.18) 

and  reduces  to  the  stress  in  a  uniaxial  specimen  when  the  right  side  of  the  equation  is  written 
for  that  condition. 

A  definition  of  equivalent  plastic  strain  is  not  quite  as  simple.  One  approach  has  been 
to  define  the  equivalent  strain  increment,  de^,  in  terms  of  an  increment  of  plastic  work,  i.e., 

dWp=aede£  (3.19) 
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Also,  we  must  have 


dWp  =  ffjj  dey  (3.20) 

which  gives 

de^  =  l/ae  CTjj  dejj  (3.21) 

It  was  stated  previously  that  hydrostatic  pressure  has  no  appreciable  effect  on  inelastic 
deformation  and  thus  that  there  is  no  change  in  volume  as  plastic  flow  occurs.  Expressed 
mathematically,  this  means  that  the  trace  of  the  strain  tensor  vanishes,  that  is, 

eH  =  0  (3.22) 

and  the  deviatoric  stress  tensor  may  be  substituted  for  the  stress  tensor  in  equation  (3.21): 

def  =  l/ae  ajjde?  (3.23) 

An  explicit  expression  for  the  equivalent  plastic  strain  may  now  be  obtained  by  taking 
the  scalar  product  of  equation  (3.23)  with  itself  obtaining 

depe  dee  =  V°l  aij  aij  defj  deij  (3-24) 

By  using  equation  (3.18),  the  expression  for  the  equivalent  plastic  strain  now  becomes 

def=  (2/3  defj  deg)1'*  (3.25) 

3.7  STRESS-TO-STRAIN  TRANSFORMATION 


With  the  above  expressions  for  equivalent  stress  and  strain,  we  may  now  look  at  the 
details  of  each  of  the  hardening  rules  that  we  will  be  considering.  In  the  isotropic-hardening 
rule,  the  change  in  the  loading  surface  is  reflected  through  a  change  in  effective  stress  which  is 
a  function  of  the  effective  strain,  that  is,  we  now  have 

f  =  J2(aij)-l/3ae(eP)2  (3.26) 

as  the  expression  for  the  loading  function. 

At  a  material  point  where  an  increment  of  plastic  strain  occurs,  equations  (3.14), 
(3.15),  and  (3.17)  must  be  true.  The  derivatives  of  the  loading  function  with  respect  to  stress 
and  plastic  strain  are  needed  in  order  to  determine  the  explicit  expression  for  equation  (3.17) 
for  isotropic  hardening.  Thus,  we  can  obtain  the  following 


9f  _  8f  dJ2_  , 

9(7:;  9J?  9(7;; 


(3.27) 
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(3.28) 


_9f  =  9£^_e^=  2j3g 
aefj  8ae  3e£  3efj  6  3e{!  3ejj 

The  derivative  of  the  loading  function  with  respect  to  plastic  strain  is  then  obtained  by  substi¬ 
tuting  equation  (3.21)  into  (3.28): 


9f-  o/,  3CTe 

=  -2/3  —  o. 


3ep. 

ij 


3e 


p  «J 


(3.29) 


Thus,  by  using  equation  (3.18),  the  following  becomes  the  linear  relationship  for  an  increment 
of  plastic  strain  in  terms  of  an  increment  of  stress: 


de 


p  _ 
mn 


amn  °k2 


3ae 
2/3  3^ 


dakC" 


pmn  °k£ 
3  a 

4/9  —  o\ 
3ep  6 


da 


kC 


(3.30) 


The  linearity  of  the  factor  on  dakg  is  now  evident  since  it  is  a  function  of  only  the  current 
deviatoric  stress,  the  current  effective  stress,  and  the  current  slope  of  the  effective  stress- 
effective  plastic  strain  curve. 

The  loading  function  for  the  kinematic-hardening  rule  must  be  different  from  equation 
(3.26)  since  the  loading  surface  is  now  allowed  to  shift  but  not  to  grow  in  size.  Incorporating 
these  two  requirements,  the  loading  function  becomes 

f  =  J2  (°ij  _aij) "  !/3  ao  (3-3D 

The  array  ay  is  the  shift  of  the  center  of  the  loading  surface  with  changing  plastic  strain, 
that  is, 

day  =CdePj  (3.32) 

where  C  is  a  constant  factor  dependent  on  the  effective  stress -effective  plastic  strain  curve.  As 
before,  aQ  is  the  yield  strength  of  the  material  which  is  obtained  in  a  uniaxial  test.  The  above 
form  of  the  kinematic  hardening  is  that  proposed  by  Shield  and  Ziegler  [56]. 

To  reemphasize  the  difference  between  equations  (3.26)  and  3.31),  the  dependence  of 
the  loading  function  on  the  plastic  strain  for  isotropic  hardening  comes  through  the  change  in 
effective  stress  which  is  a  function  of  the  effective  plastic  strain  whereas  for  kinematic  harden¬ 
ing,  the  change  in  plastic  . strain  causes  the  center  of  the  loading  surface  to  shift  by  the  amount 
specified  by  equation  (3.32). 

Equations  (3.14),  (3.15),  and  (3.17)  must  be  satisfied  for  kinematic  hardening  as  they 
were  for  the  isotropic  hardening,  Instead  of  equations  (3.27)  and  (3.28),  we  obtain  the  follow¬ 
ing  two  expressions 
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(3.33) 


_3f  _  i  t 
aajj  CTii  ~  iJ 


where 


_3fs  cii 

9eS  0CTij 


aij  aij  ^ij  akk 


Substituting  equations  (3.33)  and  (3.34)  into  equation  (3.17),  we  obtain 

(°H-a;s)(o'  n-a'  ) 


P  _  v  ij  ij/v  mn  mn'  3 

emn=CKrak£)(ffk£-akj)  °ij 


(3.34) 


(3.35) 


(3.36) 


The  denominator  in  the  above  expression  may  be  simplified  since  as  already  stated,  the  size  of 
the  loading  surface  does  not  change  and  thus  we  may  write 


(akC-akc)Ke-ake)=2/3oo 


(3.37) 


Thus,  equation  (3.36)  may  be  written  as 


p  _(qij  -alj)(omn  -amn> 
m"  2/3  Co\  " 


(3.38) 


The  problem  now  is  to  determine  the  value  of  C  in  the  above  equation.  This  can  be 
done  by  examining  the  form  of  the  equation  for  a  uniaxial  loading  test  specimen.  Equation 
(3.38)  becomes 


dell  = 


(2/3  cfj j  -a^j)  (2/3  Ojj-a  jj) 


2/3  C  a2 


don=  3ae 


(3.39) 


Written  for  this  case,  equation  (3.37)  becomes 

(2/3  (Xjj  a 1 1 ) (2/3  CTjj  -  Hjj)  +  (-1/3(7j|  -a22)(_l/3c7jj  ~a22^ 
■+■  (— 1/3  CTj  833 )  (— 1/3  <7  j  j  -833)  =  2/3  o2 

The  incompressibility  of  plastic  deformations  requires  that 

a22  =  a33  =“  1/2  a'jj 

which  can  be  substituted  into  equation  (3.40)  to  obtain 


(3.40) 


(3.41) 


3/2(2/3a11-a,11)2=2/3o2 


(3.42) 


29 


Finally,  combining  equation  (3.42)  with  equation  (3.39),  we  obtain 


C  =  2/3 


9eep 


(3.43) 


As  proposed  by  Shield  and  Ziegler  [56]  and  formulated  by  equations  (3.31)  and  (3.32), 
the  form  of  kinematic  hardening  was  intended  for  linear  hardening  where  the  slope  of  the 
equivalent  stress-equivalent  plastic  strain  curve  is  a  constant.  Nothing  in  the  above  derivation 
prohibits  modeling  the  equivalent  stress -equivalent  plastic  strain  curve  as  a  series  of  straight 
lines  where  the  slope  would  be  constant  for  each  segment.  This  idealization  is  shown  in 
Figure  3.3.  Thus,  we  may  assume  that  the  above  derivation  is  applicable  in  each  one  of  the 
straight  line  segments  of  the  idealized  curve. 

Now  substituting  equation  (3.43)  into  equation  (3.38),  we  obtain  the  desired  relation¬ 
ship  for  kinematic  hardening: 


de 


_  (gij  aij  ^°mn  amn^ 


mn 


4/9 


3ae  2 

_—CT° 

3ep 


daij 


(3.44) 


The  similarity  is  quite  interesting  when  the  above  equation  is  compared  with  the  equiva¬ 
lent  one  for  isotropic  hardening,  equation  (3.30).  The  kinematic-hardening  relationship  can  be 
obtained  from  that  for  isotropic  hardening  by  simply  substituting  the  shifted,  i.e.,  the  corrected 
deviatoric  stress  tensor,  for  the  usual  deviatoric  stress  tensor  and  using  the  intial  proportional 
limit  of  the  material  instead  of  the  current  value  of  effective  stress. 


3.8  STRAIN -TO -STRESS  TRANSFORMATION 

We  now  proceed  to  derive  the  stress -strain  relationships  which  will  be  required  for  the 
finite  element  method.  Equations  (3.30)  for  isotropic  hardening  and  equation  (3.44)  for 
kinematic  hardening  transform  an  increment  of  stress  to  an  increment  of  plastic  strain.  An 
inverse  relationship  which  transforms  total  strain  to  stress  is  required.  The  derivation  will  be 
carried  out  in  terms  of  isotropic  hardening;  the  transition  to  kinematic  hardening  is  easily  per¬ 
formed  as  described  in  the  previous  paragraph. 

The  assumption  is  first  made  that  strain  may  be  decomposed  as  follows 

de”  =  de^T  dejj  (3.45) 

T  F 

where  de^  and  de^1  are  the  total  and  elastic  strains,  respectively.  We  can  then  write 

daij=Cijk£(dek£“ak£A)  (3-46) 
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Figure  3.3  -  Idealization  of  Plastic  Stress-Strain  Curve 
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where  C^g  is  the  elastic  constitutive  tensor  which  contains  Young’s  modulus  and  Poisson’s 
ratio.  Equation  (3.13)  has  been  used  in  writing  the  above  equation. 

From  equation  (3.30),  an  explicit  expression  for  A  may  be  written: 


A  = 


ij 


a  daii 
3ae  o  J 

4/9 -V 

3*J 


(3.47) 


Now  multiplying  equation  (3.46)  by  ajj  and  combining  with  equation  (3.47),  we  may  find  an 
expression  for  A  which  is  a  function  of  strain  and  the  current  state: 


A  = 


°ij  Cijk£  dek£ 


3ffe  o 

4/9  3e^  +  aij  CijkE  akC 


(3.48) 


Equation  (3.48)  is  now  substituted  into  equation  (3.46)  to  obtain  the  desired  result: 

^ijkC  apq  Cpqrs 


daij  = 


C  •  - 
ijrs 


da 


4/9  °e  +  apq  Cpqrs  ars 


deT 

rs 


(3.49) 


The  implementation  of  the  above  equation  will  be  given  later. 

In  addition  to  equation  (3.49),  an  explicit  expression  for  the  differential  of  the  equiva¬ 
lent  plastic  strain  will  be  required  during  the  solution  of  problems.  This  can  be  written  as 


2/3  ge  P*j  C9rs  dcT 
3ae 

4/9r7ae+CTijCijk£ak£ 

3ee 


(3.50) 


► 
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CHAPTER  4 
GEOMETRIC  DETAILS 


The  development  of  an  adequate  shell  element  requires  that  considerable  attention  be 
given  to  how  the  geometry  of  the  shell  surface  is  used  in  the  formation  and  implementation  of 
the  element.  As  pointed  out  in  the  first  chapter,  Gallagher  [15]  has  indicated  that  this  has 
been  one  of  the  most  difficult  hurdles  to  overcome  in  developing  shell  elements. 

4.1  SURFACE  DEFINITION 

The  formation  of  the  stiffness  matrix  for  the  quadrilateral,  the  determination  of  the 
loading  on  the  element,  and  the  transformation  of  the  stiffness  matrix,  the  displacements,  and 
the  loads  require  detailed  information  on  shell  geometry. 

The  theory  for  the  element  has  been  developed  for  application  to  surfaces  that  can  be 
represented  by  a  quadratic  form.  Referring  to  Figure  4.1,  the  surface  must  be  able  to  be 
expressed  by  the  function  f  where  according  to  McConnell  [57], 

f(x')  =  QijX'xj  +  bkxjc  +  C  =  0  (4.1) 

Here  x'  is  the  auxiliary  coordinate  system  of  the  surface  and  and  bk  are  constants. 

Cartesian  coordinate  systems  shown  in  the  figure  are  the  local  element  coordinate  system  and 
the  global  coordinate  system.  Each  element  will  have  its  own  coordinate  system  as  will  each 
surface. 

It  will  also  be  necessary  to  represent  the  surface  by  a  vector  function 

r  =  Xj  ((3pj32)ii  +  x2  (j3 1 ,  /32 )  ^  +  x3  (/Jp^)^  (4.2) 

where  f and  ^  are  the  Gaussian  coordinates  on  the  surface,  not  necessarily  orthogonal.  The 
above  is  a  parametric  representation  of  the  surface  in  terms  of  the  coordinates  |3j  and  (3 2  ■ 

The  quadratic  form  equation  (4.1),  and  the  parametric  representation,  equation  (4.2), 
are  given  in  the  appendix  for  the  surfaces  that  are  treated  by  the  theory  developed  here. 

4.2  ELEMENT  DEFINITION  AND  BOUNDARY 

The  element  corners  are  given  by  nodal  points  1,  2,  3,  and  4  shown  in  Figure  4.2. 

The  local  element  coordinate  system,  x;,  is  established  by  the  first  three  nodes.  In  other 
words,  these  nodes  set  the  X-Y  plane  and  the  fourth  node  is  generally  not  in  the  plane.  The 
element  sides  must  be  defined  in  a  unique  way  so  that  there  will  be  no  overlapping  or  gaps 
between  adjacent  elements.  This  is  done  by  letting  the  side  of  the  element  be  along  a  line 
which  is  represented  by  a  linear  relationship  between  the  coordinates,  that  is,  each  of 

the  sides  is  given  by  the  equation 
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quadric  surface 


Figure  4.2  -  Element  in  Local  System 
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where  ctj  is  a  constant  and  i  assumes  the  values  from  one  to  four. 

The  midside  nodes  lie  on  the  above-defined  lines.  The  coordinates  of  the 

nodes  on  each  side  are  the  mean  of  the  (/Jj,/^)  coordinates  of  the  two  nodes  which  define 
the  side. 


4.3  SLOPE  OF  SURFACE  IN  ELEMENT  SYSTEM 

During  the  computation  of  the  element  stiffness  and  the  distributed  loading  on  the 
element,  the  slope  of  the  surface  is  required  in  the  local  element  system.  This  is  obtained  by 
expressing  the  quadratic  form  of  equation  (4.1)  in  terms  of  the  local  element  system.  We  use 
the  transformation  relationship  between  the  element  system  and  the  coordinate  system  of  the 
quadratic  surface.  This  is  written  as 

xi  =  aikxk  +  xi  (4-3> 

where  the  origin  of  the  element  system,  nodal  point  1 ,  is  given  by  the  vector 

R'  =  x[x|  (4.4) 


The  transformation  array,  aik,  is  determined  by  using  the  first  three  nodal  point  coordinates  as 
described  above. 

Now  substituting  (4.3)  into  (4.1),  we  obtain  an  expression  for  the  surface  in  terms  of 
the  element  system: 


1,  1, 

f(x)  =  Qy  aik  ajk  xk  x£  +  Qij  aik  xR  Xj  +  Qy  ajC  x£  Xj 


1,1, 


(4.5) 


+  Qij  xixj  +  bmamnxn  +  bmxm+C  =  0 


The  desired  slopes  are  the  derivatives  of  z  with  respect  to  x  and  with  respect  to  y. 
Since  it  would  be  extremely  tedious  to  solve  for  an  explicit  expression  for  z  in  terms  of  x  and 
y  (one  would  have  to  be  derived  for  each  surface),  we  choose  to  use  implicit  differentiation 
and  tensor  calculus  to  obtain  the  required  slopes.  As  an  example,  to  compute  the  derivative 
of  z  with  respect  to  x,  we  can  first  write  the  derivative  of  f  with  respect  to  x 


9f  +  9f  9z 
9x  9z  9x 


Thus  solving  for  the  slope,  we  obtain 


9f 

9z  9x_ 
9x  "9f 
9z 


(4.6) 


(4.7) 
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Thus  the  derivatives  of  the  function  that  expresses  the  surface  in  the  element  system  are  re¬ 
quired.  Differentiating  equation  (4.5),  we  obtain 


3f  1 

9^  Qij  aik  aj£  ^ko  xfi  +  Qj  aik  aj£  xk  +  Qij  aik  ^ko  xj 

1, 

+  Qijaje6Co  xi  +  bm  amn5no 

9f  „  1, 

9X  -  2Qij  aim  ajC  x£  +  2Qij  aim  xj  +  bp  apm 

Now  substituting  (4.8)  into  (4.7),  we  obtain  the  slope  in  the  x-direction  as 


9z 

9x 


2Qijaiiajgx£  +  2Qjjaiixj  +  bP  aPi 

2Qij  ai3ajfix£  +  2Qijai3xj  +  bp  ap3 


and,  similarly,  the  slope  in  the  y-direction  as 


9z 

9y 


2Qjj  ai2  ajg  XC  +  2Qjj  ai2xj+bp  ap2 

2Qijai3aj£xe  +  2Qij  ai3  xj  +  bpap3 


(4.8) 


(4.9) 


(4.10) 


Thus  equations  (4.9)  and  (4.10)  are  evaluated  at  each  of  the  integration  points  in  the 
element  and  then  used  to  form  the  stiffness  of  the  element  and  to  compute  the  consistent 
nodal  forces  that  are  equivalent  to  the  distributed  loading.  It  should  be  noted  that  at  a  par¬ 
ticular  integration  point  in  the  x-y  plane,  everything  is  known  that  is  required  in  the  equa¬ 
tions  except  the  value  of  z,  the  distance  of  the  element  above  or  below  the  x-y  plane.  The 
value  of  z  can  be  obtained  by  using  the  x-y  coordinates  of  the  integration  point  in  equation 
(4.5),  which  can  then  be  reduced  to  a  quadratic  equation  in  z.  The  correct  value  of  z  will  be 
the  one  which  has  the  least  absolute  value. 

Although  all  of  the  above  may  seem  to  be  somewhat  involved,  the  only  other  alterna¬ 
tive  would  be  to  assume  that  the  height  of  the  surface  above  the  x-y  plane  can  be  represented 
adequately  by  a  polynomial,  a  cubic  for  example.  Thus  since  there  are  then  ten  free  parameters 
in  a  cubic  polynomial,  ten  values  of  the  height  would  have  to  be  computed  and  then  ten 
equations  would  have  to  be  written  of  the  type 


2  2  3 

O  I  0  ,  O  ,  O “  .  00,  o^  ,  o^ 

z  =aj+a2X  +  a3y  +a4x  +  a3x  y  +a^y  +ayX 


2 

I  O  O,  00“".  O'" 

+  a8x  y  +a<,x  y  +a]0y 


(4.11) 


Here  (x°,  y°)  are  the  coordinates  of  the  point  where  the  height  z°  is  evaluated.  The  equations 
would  then  have  to  be  solved  for  the  coefficients  aj.  Finally,  the  value  of  the  slopes  at  each  of 
the  integration  points  could  then  be  calculated. 
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It  is  not  clear  that  the  cubic  fit  method  is  more  efficient  than  the  previous  devised 
method  that  is  used  in  the  quadrilateral  element.  However,  it  is  certain  that  the  previous 
method  is  more  accurate  than  the  cubic  fit  since  it  is  exact.  The  derivatives  of  functions  that 
are  represented  by  polynomial  fits  are  notoriously  poor.  Thus  one  would  not  expect  the 
slopes  obtained  from  the  cubic  approximation  to  be  very  accurate. 


4.4  SURFACE  VECTORS 


For  the  calculation  of  the  element  stiffness  and  the  transformations  required  in  imple¬ 
menting  the  element,  it  will  be  convenient  to  be  able  to  compute  vectors  on  a  line  on  the 
surface.  We  will  assume  that  the  line  can  be  represented  by  a  linear  function  of  0]  and  0-,. 
Thus  the  previously  described  element  boundary  fits  this  description.  Moreover,  we  will 
further  assume  that  the  surface  boundary  can  be  adequately  represented  by  a  series  of  these 
linear  curves.  The  convenience  of  having  the  capability  to  compute  these  vectors  will  become 
apparent  later  when  we  describe  the  computation  of  the  element  stiffness  and  the  transforma¬ 
tions. 

Consider  a  portion  of  a  surface  shown  in  Figure  4.3  which  has  been  modeled  by  a 
number  of  the  quadrilateral  finite  elements.  The  line  SI  is  a  linear  curve  between  two  ele¬ 
ments  and  the  line  S2  is  a  portion  of  the  boundary.  Both  lines  are  linear  in  the  sense 
described  above. 

The  vectors  b,  t,  and  n  are  unit  vectors  at  the  particular  points  J  and  L.  The  vector  n 
is  normal  to  the  surface  and  is  computed  by  first  taking  the  cross  product 


dr  dr 
—  x  =  n 
9/3j  002  ~ 


(4.12) 


where  r  is  the  vector  function  that  defines  the  surface  and  is  given  in  equation  (4.2).  The  unit 
normal  is  then  obtained  by  dividing  the  normal  vector  by  its  magnitude,  i.e., 


(4.13) 


To  compute  the  tangent  vector,  we  must  first  express  the  line  as  a  vector  function.  We 
do  this  by  using  the  equation  of  the  line  which  is  a  linear  function  of  0j  and  02.  Thus  solv¬ 
ing,  say,  for  02  in  terms  of  0j  and  substituting  into  equation  (4.2),  we  obtain 

r(0j ,  02 (0j ) )  =  P(0j )  =  F(0j  )i'j+  GOSj  )i2  +  H(0j  )i^  (4. 1 4) 


where  p  is  now  the  vector  equation  of  the  line  between  the  elements  or  of  the  segment  of  the 
boundary.  This  depends  on  the  particular  situation  or  case  for  which  the  computation  is  being 
performed.  We  now  obtain  the  tangent  vector  to  the  line  by  computing  the  derivative  of  p  with 
respect  to  0  j  to  obtain 
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Figure  4.3  —  Surface  Vectors 


(4.15) 


9 p  3F  A  dG.  3Ha 

wr*/'**^**^ 

The  unit  tangent  vector  is  now  obtained  by  dividing  the  above  vector  by  its  length,  as  was 
done  for  the  normal  vector.  The  positive  sense  of  the  tangent  vector  is  taken  in  the  direction 
of  increasing  j3j.  If  the  line  happens  to  be  on  a  constant  line,  the  above  process  for  com¬ 
puting  the  tangent  vector  is  performed  by  using  ^  as  the  parameter. 

The  final  vector  b  is  obtained  by  performing  the  cross  product 

b  =  txn  (4.16) 

This  vector  is  a  unit  vector  since  the  two  vectors  used  in  the  cross  product  are  unit  vectors. 

a 

However,  the  b  vector  is  not  necessarily  the  binormal  to  the  curve  between  two  elements  or 
the  curve  along  the  boundary,  but  it  is  a  unique  vector  which  is  sufficient  for  our  purpose, 
namely,  to  compute  the  element  stiffness  and  perform  transformations.  If  it  so  happened 
that  the  linear  line  on  the  surface  was  also  a  geodesic  line,  then  the  normal  to  the  curve  would 
coincide  with  the  normal  to  the  surface  and  then  the  vector  would  be  the  binormal  to  the 
curve  at  the  point  in  question. 
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CHAPTER  5 

QUADRILATERAL  SHELL  FINITE  ELEMENT 

The  derivation  of  the  shell  quadrilateral  finite  element  will  be  obtained  by  showing  the 
particular  form  of  the  basic  formulation  presented  in  Chapter  2.  At  particular  points  in  the 
derivation,  we  will  draw  on  the  geometric  information  given  in  Chapter  4. 

5.1  STRAIN  DISPLACEMENT  RELATIONSHIP 

The  strain  displacement  relationship  for  the  quadrilateral  element  is  that  due  to 
Marguerre  [58]  and  also  presented  in  Washizu  [38].  It  is  assumed  that  the  strain  is  repre¬ 
sented  adequately  by  letting  the  element  be  a  segment  of  a  shell  that  does  not  differ  greatly 
from  a  segment  of  a  flat  plate  which  has  an  initial  deformation.  This  approximation  will  not 
limit  the  application  of  the  element  and  falls  within  the  well-known  shallow  shell  theory. 

Thus  the  only  limitation  imposed  is  that  each  individual  element  is  a  shallow  shell.  The 
specifics  and  source  of  this  limitation  will  be  brought  out  in  the  derivation. 

Referring  to  Figure  5.1.,  the  undeformed  shell  element  is  assumed  to  be  represented 
by  the  mapping 

x  =  X 

y  =  Y  (5.1) 

z  =  w 

where  w  is  the  distance  from  the  (X,  Y)  plane  to  the  midfiber  of  the  shell.  After  load- 
induced  deformation,  the  midfiber  mapping  becomes 

x  =  X  +  u 

y  =  Y  +  v  (5.2) 

z-  w  +  w 

where  we  consider  u,  v,  w,  and  w  to  be  functions  of  X  and  Y  only. 

Prior  to  load -induced  deformation,  a  point  on  the  midfiber  can  be  represented  by  the 

vector 

r^0)  =  Xij  +  Yi2  +  wi3  (5.3) 

and  a  point  distance  £  above  or  below  the  midfiber  by 

£(o)=I(o°)+£n(o)  (5.5) 

where  n(0^  is  the  normal  to  the  surface  before  deformation.  This  normal  is  computed  by 
taking  the  cross  product 
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Figure  5.1  -  Mapping  of  Shell  Surface 
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3r(0)  ar(o) 

~o  ~o 

-  X  - 


fj(o)  =  9X  8Y  (5 

ar(0)  3r(o) 

-a-  x  -S- 

ax  3Y 

After  load  is  applied  to  the  element,  a  point  on  the  midsurface  assumes  the  position 


rQ  =(X  +  u)ij  +(Y  +  v)i2  +(w  +  w)i3  (5.6) 

The  normal  to  the  new  position  is  now  computed  by 

a_Epyaio 

ax  aY 

The  approximation  of  the  shallow  shell  theory  in  Cartesian  coordinates  is  now  intro¬ 
duced  through  the  assumption  that 

(i)2.(nj.im«>  <5.8, 

and  thus  sets  the  limitations  on  how  “deep”  the  shell  element  can  be  without  degrading  the 
results  of  the  computations. 

We  will  further  assume  that  the  shell  element  is  thin  and  the  curvature  is  not  very  large, 
such  that 


u  a2w  y  a2w  ..  a2w 
ax2  ’  3Y2  ’  axaY 


« 1 


and  that  multiples  of  £  cause  lower  order  magnitude  terms.  Also,  since  we  will  be  using  the 
approximation  that  the  stretching  terms  are  small,  we  may  assume  that 


3u  8u  3v_  3v  . 

ax’aY’ax’aY’ 


(5.10) 


Imposing  (5.8),  the  normal  to  the  element  surface  before  loading  may  be  computed  by 
(5.5)  to  obtain 


„(o)  _3w o  * 

ax  1  9Y  2  ’3 


(5.11) 


After  deformation,  the  normal  to  the  midsurface  is  obtained  from  (5.6) 

a  /3w  ,  3wV  /3w  ,  aw\-  ,  c 

n '"'lax  axh  "VaY+aYj12+13 


(5.12) 
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where  (5.10)  has  been  used  along  with  the  assumption  that  the  square  of  the  derivative  of  w 
may  be  ignored  in  computing  the  unit  vector. 

The  additional  constraint  of  the  Kirchoff  hypothesis  will  be  imposed  by  the  equation 

r  =  r0  +  £n  (5.13) 

which  is  the  position  vector  of  any  point  through  the  thickness  after  deformation. 

We  may  now  compute  the  strain  tensor  by  using  (2.3).  Before  proceeding,  it  is  impor¬ 
tant  to  note  that  equation  (2.3)  is  based  on  the  assumption  that  the  line  segment  before 
deformation  can  be  written  directly  in  terms  of  the  increment  of  the  Lagrangian  coordinates. 
This  is  not  true  for  the  case  here  where  the  actual  load -induced  deformation  starts  out  from  a 
deformed  configuration,  i.e.,  the  element  shell  surface.  We  may  still  use  a  Lagrangian  descrip¬ 
tion  as  long  as  all  quantities  are  referenced  to  the  basic  configuration  in  the  (X,  Y)  plane. 

Thus,  using  (5.4)  and  (5.13),  we  may  write  (2.3)  in  the  form 


.  .  /3r  \/3r  \  /3r<0>  \  /3r<°>  ' 

ds  -dS  -^9X.dXij-\axjdXjj“\3Xi  dXi)*\3Xj  dXi 


(5.14) 


=  2EijdXidXj 


The  final  expressions  for  the  strain  quantities  are  now  obtained  by  substituting  the  previous 
expressions  into  (5.14)  and  dropping  all  terms  which  are  assumed  to  be  small  for  the  reasons 
stated  previously.  For  any  position  through  the  thickness  of  the  shell,  the  three  strain  terms 
become 


F  _  du  3w  3w  .  /3wV  _  y  32w 
fcxx  ax  3X3X  1  \dXj  50X2 

f  3v  3w  3w  .  ,?/3wV  t  w 
tYY=3Y  3Y3Y  1/Z\3Y/  ~%Y2 


(5.15) 


(5.16) 


iv  _  3u  3v  3w  3w  3w  3w  3w  3w  _  ?  32 w 
ztXY  3Y  3X  3X3Y  3Y3X  3X3Y  5  3X3Y 


(5.17) 


The  above  equations  are  appropriate  for  small  direct  strains  and  moderate  rotations. 


5.2  ASSUMED  DISPLACEMENTS 

The  quadrilateral  element  on  the  surface  of  the  shell  was  shown  in  Figure  4.2.  We 
must  now  examine  the  region  that  is  the  projection  of  the  element  onto  the  local  (X,  Y) 
plane.  In  the  Lagrangian  sense,  this  region  contains  the  material  coordinates  in  which  we  will 
assume  that  a  certain  displacement  exists  and  over  which  we  will  have  to  express  the  loading 
and  perform  the  integration  for  the  stiffness  equation. 
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On  the  shell  surface,  the  element  is  quadrilateral  in  shape  but  the  projected  region  in 
the  (X,  Y)  plane  is  not.  Nodal  points  1  -3  establish  the  (X,  Y)  plane  and  thus  the  local 
element  coordinate  system.  The  remaining  nodes  are  projected  onto  the  plane  to  establish  the 
eight-sided  polygon  shown  in  Figure  5.2.  The  true  region  is  curvilinear  in  shape  for  the 
general  application,  but  the  above  approach  should  be  sufficient  to  represent  the  area. 

Following  the  work  of  Fraeijs  de  Veubeke  [59],  we  use  patching-type  bicubic  functions 
for  the  assumed  displacement  pattern.  This  approach  was  used  originally  for  studying  plate 
problems,  but  one  of  the  requirements  to  maintain  compatibility  between  shell  elements  is 
that  the  same  order  (here  bicubic)  polynomial  must  be  used  for  each  of  the  three  u,  v,  and  w 
displacements.  We  refer  here  only  to  the  displacements  of  the  midfiber  of  the  shell  element 
and  consider  only  the  u  displacement  first. 

Referring  to  Figure  5.2,  the  complete  ten-term  cubic  polynomial  is  used  over  the 
entire  region,  that  is, 


uc=3j  +a2X  +  a3Y  +  a4X2+a5XY  +  a6Y2+a7X3  +  a8X2Y+a9XY2  +  a10Y3  (5.18) 

Next,  the  two  auxiliary  coordinate  systems  shown  in  the  figure  are  introduced.  An  additional 
function 


u'=an  Y'2+a12Y'3+a13X'Y'2  (5.19) 

is  defined  to  be  acting  in  the  region  where  Y'  is  greater  than  zero.  Also,  a  third  function 

u"  =a14  Y"2  +  a15  Y"3  +a16X”  Y"2  (5.20) 


is  defined  to  be  acting  in  the  region  where  Y”  is  greater  than  zero.  Each  of  the  functions  is 
written  with  unknown  coefficients  which  are  to  be  determined.  Thus,  the  u  displacement  at 
any  point  within  the  region  can  be  expressed  by 


u  =  uc  +  u'  +  u”  (5.21) 

with  the  prime  and  double-prime  functions  being  defined  as  above.  All  three  displacements  on 
the  right  side  of  (5.21)  are  directed  parallel  to  the  positive  X  axis. 

Similar  to  the  u  displacement,  the  v  and  w  displacements  are  defined  by 


and 


v  =  vc  +  v'  +  v'' 


(5.22) 


w  =  wc  +  w'  +  w"  (5.23) 

where  the  prime  and  double-prime  terms  have  the  same  meaning,  only  now  the  range  on  the 
coefficient  subscripts  is  17  to  32  on  v  and  33  to  48  on  w.  We  may  express  this  in  matrix 
form,  which  is  the  equivalent  of  equation  (2.18): 
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Figure  5.2  —  Auxiliary  Coordinate  Systems 
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(5.24) 


where 


|r]  =[1  X  Y  X2  XY  Y2  X2  X2Y  XY2  Y3  Y'2  Y'3  XY‘  Y"J  X"Y 


//2  \r  f  \y  1 1  5 


rt\ 2 


(5.25) 


The  above  displacements  are  only  for  the  midfiber  of  the  shell  element  where  equation  (2.18) 
was  written  for  the  general  case  for  any  point  in  the  body.  The  displacement  for  any  point 
through  the  thickness  of  the  element  will  come  from  the  constraints  of  the  Kirchoff  hypothesis, 
equation  (5.13). 


5.3  GENERALIZED  DISPLACEMENTS 


The  concept  of  generalized  displacements  was  introduced  in  Chapter  2  through  equation 
(2.20)  and  the  characteristics  of  these  displacements  were  stated.  One  characteristic  was  that 
the  number  of  generalized  displacements  is  equal  to  the  number  of  undetermined  coefficients. 
Thus  for  the  quadrilateral,  we  will  introduce  48  generalized  displacements  in  accordance  with 
equation  (5.24).  Our  eventual  goal  here  is  to  express  the  displacements  of  the  midsurface  of 
the  element,  equation  (5.24),  in  terms  of  the  generalized  displacements. 

The  nine  generalized  displacements 


IvJ 


9u  9u_  9v  9v_  9w9w 
U  V  W  9X  9Y  9X  9Y  9X  9Y_ 


(5.26) 


will  be  used  at  nodes  1,  2,  3,  and  4,  and  the  three  generalized  displacements 


(5.27) 


will  be  used  at  each  of  the  midside  nodes  5,  6,  7,  and  8,  giving  the  required  total  of  48.  The 
b  direction  in  equation  (5.27)  is  the  direction  defined  by  equation  (4.16)  and  shown  at  point 
J  in  Figure  4.3. 

It  was  pointed  out  in  Chapter  4  that  the  direction  and  positive  sense  of  b  in  equation 
(4. 1 6)  is  unique  and  is  referenced  to  the  coordinate  system  on  the  shell  surface.  Thus  the 
derivative  of  the  displacements  with  respect  to  b  in  two  adjacent  elements  will  have  a  positive 
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sense  in  the  outward  direction  in  one  and  an  inward  positive  sense  in  the  other.  This  is  im¬ 
portant  for  maintaining  compatibility  between  the  two  adjacent  elements  since  the  displace¬ 
ments  along  their  common  boundary  line  are  expressed  in  terms  of  the  generalized  displace¬ 
ments  at  the  nodes  on  the  line. 

Expressions  for  the  generalized  displacements,  equation  (5.26)  and  (5.27),  can  be 
determined  by  using  equations  (5.24)  and  (5.25).  The  first  three  displacements  of  (5.26)  are 
the  same  as  (5.24).  The  remaining  generalized  displacements  can  be  determined  by  the  chain 
rule.  As  an  example,  from  (5.21),  we  have 


3u  3  3u^  9u^ 

ax;  aXj+  aJ>  ax:+cjiaxj' 


(5.28) 


where 


ji  3X; 


(5.29) 


and 


C:i 


dX” 
'J1  =  a3T 


(5.30) 


and  thus  the  derivatives  of  the  total  u  function  with  respect  to  X  and  Y  can  be  determined. 
The  second  order  tensors  of  (5.29)  and  (5.30)  are  arrays  of  direction  cosines  which  relate  the 
two  additional  coordinate  systems  to  the  basic  element  coordinate  system.  Also,  the  deriva¬ 
tive  of  u  at  midside  node  m  with  respect  to  the  direction  b  can  be  determined  from 


where 


3u  9uc 

~  =df3+d  k 

3bm  3Xi 


a  —  +  dm  C 

jk  axj  dk  cJk  ax” 


(5.31) 


ax* 

dfn  = -  (5.32) 

3bm 


Expression  (5.32)  is  a  shorthand  notation  for  the  two  direction  cosines  which  relate  the  X  and 
Y  axes  to  direction  b  at  midside  node  m. 

We  may  now  write  the  48  equations  in  the  generalized  displacements  so  that  each  is 
evaluated  at  its  respective  nodal  point  in  accordance  with  equation  (2.20),  but  since  the  same 
polynomial  form  is  used  for  each  of  the  u,  v,  and  w  displacements,  only  16  equations  have  to 
be  formed.  Thus,  considering  only  the  generalized  displacements  in  u,  we  may  write 


Ivj  = 


•-u4  | 


3u 

ax, 


3u  3u 
~"9X4  3Yj 

=  [6j  e  1 7] 


3u  j  jfu  3u 

3Y4  1  Sb5  3bg 


(5.33) 
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where  the  numerals  indicate  nodal  point  numbers.  The  Greek  alphabet  symbols  are  intro¬ 
duced  for  future  reference.  Now,  writing  (2.20)  for  the  u  displacement 


and  inverting,  we  obtain  equation  (2.21)  for  the  quadrilateral  element 


(5.34) 


(5.35) 


where  the  bar  has  been  written  over  the  A  matrix  to  indicate  that  only  the  u  displacement  has 
been  considered. 

Equation  (5.35)  expresses  the  unknown  coefficients  of  the  u  displacement  within  the 
element  in  terms  of  the  u  generalized  displacements  at  the  corner  and  midside  nodes.  The  A 
bar  matrix  is  identical  for  the  v  and  w  displacements  and  does  not  have  to  be  computed. 

The  midside  nodes  in  (5.35)  are  undesirable  and  are  difficult  to  implement  in  a  finite 
element  computer  program.  In  the  present  application,  they  are  eliminated  by  assuming  a 
linear  variation  of  u,  v,  and  w  along  the  side  of  the  element  in  the  b  direction.  Thus,  consid¬ 
ering  only  the  u  displacement,  the  derivative  of  u  in  the  b  direction  at  the  midside  node  is 
assumed  to  be  the  mean  of  the  combination  of  the  derivative  of  u  with  respect  to  X  and  Y  at 
each  of  the  nodes  which  define  the  side  of  the  element.  This  may  be  shown  by  first  rewriting 

(5.35)  in  the  form 

[|S|5  i  l51'  i  [S17]  jz  >  (5.36, 

The  linear  variation  may  be  expressed  by 

|tHg1M  (5.37) 

where  the  terms  of  G  are  factors  of  one-half  and  direction  cosines.  Thus  as  an  example,  using 
the  notation  of  (5.32),  we  may  write 


3u 

ab* 


=  1/2  d 


3u_ 

l]  ax. 


+  d' 


in  +de  in+d6  in 

1  ax3  d2  3Y2  +  d2aY3 


(5.38) 
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The  other  three  equations  are  obtained  similarly.  Substituting  (5.37)  into  (5.36), 


!  [A]e  +  [A]7  [G] 


(5.39) 


we  obtain  an  equation  for  the  coefficients  in  terms  of  the  generalized  displacements  at  only 
the  corner  nodes. 

Equation  (2.22)  may  now  be  formed  for  all  the  coefficients  and  the  generalized  dis¬ 
placements  at  each  of  the  four  corner  nodes  by  using  the  transformation  matrix  formed  in 
equation  (5.39).  Before  doing  this,  it  is  necessary  to  introduce  additional  generalized  displace¬ 
ments  which  are  the  derivatives  of  u,  v,  and  w  with  respect  to  the  Z  coordinate.  These 
additional  generalized  displacements  do  not  describe  the  deformation  of  the  element,  but  they 
will  facilitate  the  transformations  in  going  to  and  from  various  coordinate  systems  that  will 
have  to  be  performed  later.  Thus,  the  set  of  generalized  displacements  at  each  of  the  four 
corner  nodes  is  now 


Ivj 


i  =  [» 


V  w 


9u  9u  9u  9v  9v  3v  9w  9w  3w 


J, 


3X  9Y9Z  9X3Y  9Z  9X  3Y  9Z. 

The  equivalent  of  equation  (2.21)  for  the  quadrilateral  element  can  now  be  written  as 


(5.40) 


r  ^  ^ 

'v,; 

[  / 

i 

>  -  IA1  < 

v-  i 

: 

v3 

00 

03 

— 

l?4  J 

(5.41) 


where  the  columns  of  the  A  matrix  are  obtained  from  equation  (5.39)  and  are  ordered  in  accord¬ 
ance  with  the  ordering  of  the  generalized  displacement  in  (5.40).  The  columns  which  correspond 
to  the  derivatives  with  respect  to  Z  are  void. 

The  final  equation  which  expresses  the  displacements  within  the  element  in  terms  of  the 
generalized  displacements  at  the  four  corner  nodes  may  now  be  written  by  substituting  (5.41) 
into  (5.24)  to  obtain 


(5.42) 
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5.4  GENERALIZED  STRESSES  AND  STRAINS 

The  stiffness  equation,  equation  (2.30),  was  derived  for  a  general  three-dimensional 
problem  and  must  be  specialized  for  plate  and  shell  problems  which  are  solved  by  using  the 
Kirchot'f  hypothesis,  equation  (5.13).  We  begin  by  introducing  the  concept  of  generalized 
stresses  and  strains. 

Equation  (3.49)  gives  the  constitutive  relationship  for  an  increment  of  stress  in  terms  of 
an  increment  of  total  strain.  Introducing  the  notation  of  Marcal  [60] ,  we  may  write 

|AS|  =  [p-]  {AE(  (5.43) 

where  the  stresses  and  strains  have  been  put  into  vector  form  and  the  matrix  p-  is  the  term  in 
brackets  in  equation  (3.49). 

The  generalized  stresses  are  introduced  through  the  equation 


2 


where  the  stresses  are  integrated  through  the  thickness,  h,  of  the  element.  The  generalized 
stresses  are  now  the  thrust  and  moment  per  unit  inch. 

The  expression  for  generalized  strains  can  be  obtained  by  first  writing  an  increment  of 
total  strain  as 

)ae|  =  (ae(  +  Z{A«c[  (5.45) 

where  the  first  term  is  the  strain  at  the  midfiber  and  the  second  term  is  the  additional  strain 
due  to  bending.  Using  (2.15)  and  (5.15)  through  (5.17),  the  two  components  of  strain  for  the 
quadrilateral  element  become 
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and 


r  *\ 

A*XX 


\ak}=( 


A  K- 


-  A 


YY 


- 


2AKxy 

v.  J 


-2A 


82w 

8X2 

82w 

8Y2 

82w 


8X8Y^ 


(5.47) 


The  engineering  definition  of  shear  strain  has  been  used  in  (5.46)  and  (5.47).  Substituting 
(5.43)  into  (5.44)  and  using  (5.45),  we  obtain  the  relationship  between  the  components  of 
generalized  strains  and  generalized  stresses: 


)  an  )  -  r  -  -  -  J  z--~~  1  f  aE  \ 

(AM)  Jh  Lz[p-]  |Z2[p~]J  (A*) 

"2 


(5.48) 


5.5  STIFFNESS  MATRIX  AND  EQUILIBRIUM  CORRECTION  VECTOR 

The  previous  definitions  of  generalized  stresses  and  generalized  strains  set  the  form  of 
the  stiffness  matrix  and  the  equilibrium  correction  vector  of  equation  (2.30).  The  D  matrix 
of  that  equation  is  the  matrix  under  the  integral  in  (5.48).  The  B1  and  B2  matrices  of  (2.30) 
are  defined  in  Table  5.1.  The  B1  matrix  is  the  increment  in  generalized  strain  and  is  thus 
identical  to  the  combination  of  (5.46)  and  (5.47).  The  B2  matrix  comes  from  the  variation 
on  the  increment  of  strain.  The  specific  evaluation  of  the  two  matrices  required  in  the 
numerical  integration  is  performed  by  using  the  expressions  for  the  u,  v,  and  w  displacements 
given  in  (5.24).  The  numerical  integration  procedure  is  described  in  the  next  chapter. 

5.6  LOADING 

The  distributed  loading  for  the  quadrilateral  element  depends  on  the  relationship  of 
the  differential  area  on  the  shell  surface  to  the  projected  differential  area  in  the  (X,  Y)  plane. 
The  two  differential  areas  are  shown  in  Figure  5.3.  An  increment  of  distributed  load 
expressed  in  terms  of  area  on  the  shell  surface  must  be  corrected  to  the  equivalent  distributed 
load  in  the  base  plane  since  this  is  where  the  integration  will  be  performed.  Moreover,  the 
change  in  the  area  on  the  shell  due  to  a  change  in  deformation  must  also  be  determined. 

The  vector  differential  area  on  the  shell  surface  can  be  expressed  as  the  cross  product 
of  the  differential  vectors  on  the  surface 


da  =  dr  x  5r 


(5.49) 
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TABLE  5.1 

B1  AND  B2  MATRICES 


[BI]  = 


3w  . 9w\  9 
+  9X  +  gA  9X/ 9X  r 


3w  3w  ,  .  3w\  3  ,  , 

9Y+3Y  gA9Y/9Y  lrJ 


9  ,  ,  I  9  ,  ,  |  /9w  3w\  9  /9w  9w\  9  ,  ,  ,  A 3w  3  .  . 

9Y  r  I  9X  r  I  \9X  +  9X /9Y  r  \9Y  +  9Y /9Y  f  +  A9X9YlfJ 


3X2 


29X3YlrJ 


1=1  g  =  1/2 

1  =  2  g=  1 


h  =  2 


[r]  is  defined  in  equation  (5.25) 


or  by  the  tensor  equivalent 


dai  =  eijk  dxj  5xk  =eijk  xj,m  dXm  xk,n5Xn 


(5.50) 


Thus,  an  expression  for  the  shell  surface  area  can  be  obtained  in  terms  of  the  base  plane  area 
by  using  (5.1)  in  (5.50),  i.e., 


,  /  9w«  9wa  ^  \ 

d'  =  \  9X1!  -9y'2  +  13/  A 


(5.51) 


which  will  be  used  shortly  to  obtain  the  consistent  loading  vector.  As  expected,  equation 
(5.51)  is  similar  to  the  equation  of  the  normal  on  the  shell  surface,  (5.1 1). 

We  now  turn  to  the  problem  of  determining  the  change  in  the  area  due  to  the  change 
in  deformation  that  occurs  in  going  from  the  intermediate  configuration  to  the  current  con¬ 
figuration.  We  may  write  the  incremental  equivalent  of  (5.2)  as 


xi  +  Axj  =  Xj  +  Uj  +  Auj 


(5.52) 


where  it  is  understood  that  the  initial  warping,  w,  is  included  in  the  component  Uj  on  the  right 
side  for  i  equal  to  three.  The  change  in  the  differential  area  can  now  be  expressed  by 
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Figure  5.3  -  Differential  Area  on  Shell  Surface  and  in  Base  Plane 


Adai  =  eijk  d(Xj  +  Axp  5(xk  +  Axk)  -  ejjk  dXj  5xk 
Substituting  (5.52)  into  (5.53),  the  change  in  area  can  be  shown  to  be 


(5.53) 


Adai  =  eijk  <5jm 


Auk ,n  +  uj,m  Auk ,  n  +  Auj  ,m  6k  n  +  Auj ,m  uk ,n 
+  Auj,m  Auk,n)dXmSXn 


(5.54) 


Now  imposing  the  assumptions  of  (5.8)  and  (5.10),  the  final  expression  for  the  change  in  area 
turns  out  to  be 


/  3w„  9w«  „  \ 

Ad~  \  Aax‘1  _A av^  +  0i3/)dA 

Ada  =  Amj  ij  d A 


(5.55) 


It  is  interesting  to  note  that  (5.55)  is  similar  to  the  difference  in  equations  (5.12)  and  (5.1 1), 
the  difference  between  the  surface  normals  after  and  before  deformation.  The  m  components 
are  the  estimated  displacements  within  the  element  and  their  computation  will  be  described 
subsequently. 

The  loading  terms,  (2.27)  and  (2.28),  may  be  determined  for  the  quadrilateral  element 
by  using  equations  (5.51)  and  (5.55).  The  relative  first-order  amplitude  of  the  two  differential 
areas  before  loading  is 


(5.56) 


The  body  force  per  unit  volume  of  the  shell  surface  is  shown  in  Figure  5.3,  which  is 
expressed  in  terms  of  the  local  element  system  as 


g  =  gjij  (5.57) 

The  incremental  part  of  (2.27)  may  now  be  written  in  matrix  form  using  (5.56)  and  the 
notation  of  (2.32) 


|AP|  =h[A]y [r]T  |Ag|CdA  (5  5g) 

A 

where  h  is  the  thickness  of  the  element  and  the  scalar  symbol  A  indicates  integration  over  the 
(X,  Y)  plane.  The  total  load  P  of  (2.32)  at  the  intermediate  configuration  is  obtained  by 
simply  adding  all  the  increments  that  have  been  computed  previously. 

The  case  of  uniform  normal  pressure  will  now  be  considered  for  the  distributed  sur¬ 
face  loading.  The  surface  traction  acting  on  the  shell  is  given  by  the  equation 

f  =  -pda  (5.59) 
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and  the  traction  can  be  expressed  in  terms  of  unit  area  in  the  (X,  Y)  plane  by  using  (5.51 ) 
and  (5.56)  to  obtain 


Fj  =  PgjC(X,  Y)dA  (5.60) 

With  (5.60),  the  incremental  term  of  (2.27)  due  to  surface  loading  is  now 

|AF}  =  -Ap[A]yVlT{e}C(X,  Y)dA  (5  6I) 

A 

where  the  components  of  the  surface  differential  area  have  been  put  into  column  vector  form. 

The  change  in  the  surface  traction  due  to  the  change  in  deformation  may  be  written 
directly  now  as 

(aFg|  =-p[A]^* [r]T  { Am  |  C(X,  Y)dA  ($  ^ 

A 


where  the  symbol  m  is  defined  in  (5.55). 

The  £  components  in  (5.61)  are  dependent  on  the  geometry  of  the  shell  whereas  from 
(5.55),  the  m  components  are  dependent  on  the  increment  of  displacements  that  occur  within 
the  element.  Using  (5.42),  we  may  write  the  m  components  as 

Anij  =  1 0 1 0 1  r  J,j  [A]  {Av}  (5.63) 

where  the  comma  is  used  to  indicate  partial  differentiation  with  respect  to  the  index  i.  The  m 
components  can  now  be  put  into  vector  form  by  writing 

{Am|  =  [q]  [A]  |Av}  (5.64) 

where  the  last  row  of  the  matrix  q  is  void  in  accordance  with  (5.55).  Now  substituting  (5.64) 
into  (5.62),  we  obtain  the  final  form  of  the  increment  of  load  due  to  the  change  in  deforma¬ 
tion  as 


{AFG}=-p[A]y*[r]TC(X,  Y)[q]dA(A]{Av(  (5.65) 

A 

In  forming  the  above  vector,  it  was  noted  previously  that  we  need  the  increment  of  displace¬ 
ment  that  occurs  in  going  from  the  intermediate  configuration  to  the  current  configuration. 
Since  these  are  the  fundamental  unknowns,  the  same  procedure  applies  here  as  was  used  for 
estimating  the  displacements  described  for  equation  (2.30). 

Similar  to  the  total  body  force  loading,  the  total  distributed  loading,  F  in  (2.32),  is 
obtained  by  adding  all  of  the  contributions  computed  and  using  (5.61)  and  (5.65)  for  all 
previous  increments. 

The  loading  terms  derived  above  are  the  generalized  forces  and  are  conjugate  to  the 
generalized  displacements  in  the  sense  that  the  product  of  the  two  produces  work  which  is  a 


55 


scalar.  Thus,  the  choice  of  a  particular  generalized  displacements  used  and  the  application  of 
the  virtual  work  theorem  produced  the  generalized  forces. 


5.7  TRANSFORMATIONS 

It  is  important  to  develop  a  consistent  and  valid  approach  for  transforming  the 
generalized  displacements  and  forces  to  and  from  various  coordinate  systems.  The  effective¬ 
ness  of  the  element  would  be  greatly  reduced  if  its  use  were  restricted  to  one  system. 

Consider  first  the  transformation  of  the  generalized  displacements  at  a  node  given  by 
(5.40).  It  is  necessary  to  obtain  a  method  that  transforms  the  displacements  between  the 
global  system  and  the  local  element  system. 

The  transformation  of  global  coordinates  (barred)  to  local  element  coordinates  will  be 
taken  as 

X^ayXj  (5.66) 

where 


aij=cos(XijXj) 

The  generalized  displacements  at  nodal  point  i  are  partitioned  as  follows 

Lu  vw!axaYazax3Y  azaxay 
=  l  V1  V2J; 


(5.67) 


(5.68) 


where,  as  before,  the  bar  over  the  displacements  indicates  the  values  at  the  nodes.  The  transfor¬ 
mation  of  the  first  type  (superscript  1)  is  straightforward  since  they  are  vector  components. 
Thus 


v1  =  a  -  V* 
i  ij  J 


(5.69) 


where  the  capital  letters  indicate  the  generalized  displacements  in  the  global  system. 

To  see  how  the  second  type  of  generalized  displacement  (superscript  2)  will  transform, 
consider  an  increment  of  the  first  type  of  generalized  displacement  in  the  local  element 
coordinate  system,  i.e., 


dv:1  =  — r  dX; 


1  aXj 


(5.70) 


Using  (5.69),  we  obtain 


d?i  =a4(aikvk)dxj 


(5.71) 
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and  expanding  the  above 


9x 


dvi  3Xo(a‘k  V^}  9Xj  dXj 


9Xg  lk  k  9Xj 


(5.72) 

Since  the  transformation  is  independent  of  the  position  in  the  coordinate  system,  we  may  write 

(5.73) 


-i  8vk 

d’'a'kaj«ax;dX) 


where  the  inverse  of  (5.66)  has  been  used.  Comparing  (5.70)  and  (5.73),  we  see  that 

8_1  .  9  V’ 

9Xj  a‘k  ai£  9X£ 

Thus  the  transformation  of  the  second  type  of  generalized  displacement  is 

.  ■  4 
•  1  * 

7ij=aikaj£^g 


(5.74) 


(5.75) 


which  is  the  transformation  for  second -order  tensors. 

The  transformation  of  the  first  type  of  generalized  force  follows  from  (5.69)  and  the 
equality  of  work  performed  in  the  local  element  system  and  the  global  system  [61].  Equating 
the  work  performed  in  the  local  system  to  that  in  the  global  system, 

f iv/=Fi1V/  (5.76) 

and  substituting  (5.69) 

f/  (ay  V 1 )  =  Fj  V*  =FjI  Vj1  (5.77) 

we  obtain 

fNj-Fj1  <5-78> 

Now  multiplying  by  the  transformation  tensor 

<i  “kj  »ij  =  *i'  Ski  =  “kj  Fj  <s'79> 

we  obtain  the  transformation  of  the  first  type  of  generalized  forces  as 

f‘=akjF)  (5.80) 

Since  the  second  type  of  generalized  displacements  transform  as  second -order  tensors, 
the  generalized  forces  must  also  be  elements  of  a  second -order  tensor  in  order  to  be  conjugate 
to  the  displacements  in  the  sense  that  the  product  of  the  generalized  displacements  and 
generalized  forces  produce  work  which  is  a  scalar.  Thus  as  was  done  for  the  first  type  of 
forces  and  displacements  in  (5.76),  we  may  now  equate  the  work  in  the  local  and  global 
coordinate  systems: 
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(5.81) 


and  the  derivation  follows  similarly  to  that  given  for  the  first  type.  Substituting  (5.75)  into 
(5.81). 


f2;  a, 


•a.  \p-  =  F2  w2  -  p/  \j l 

lj  “ini  jn  mn  lj  ij  mn  mn 


(5.82) 


and  thus 


f2  a  ^  -  F2 

ij  im  jn  1  mn 

Now  multiplying  by  the  transformation  tensors 

*ij  a£m  aim  apn  ajn  -  *ij  ^pj  -aCm  apn  ^mn 

we  finally  obtain 


(5.83) 


(5.84) 


^p  a£m  apn  ^mn 

All  of  the  above  may  now  be  put  into  matrix  notation, 
expressed  as 


(5.85) 

Equation  (5.69)  may  be 


jvH  =  [t1  ]  jv1} 


(5.86) 


where  the  notation  is  self-explanatory. 

The  second  type  of  displacement  transformation  could  be  written  as 

|v2}  =[t2]  [V2]  [t2]T  (5.87) 

but  would  cause  problems  in  the  implementation  of  the  element.  The  element  has  been 
exercised  in  a  program  (to  be  described  in  the  next  chapter)  which  presently  has  provision 
only  for  the  transformation  of  generalized  displacements  and  forces  which  transform  as 
vectors.  The  reason  is  that  there  has  been  no  need  for  any  other  type  of  transformation  since 
the  degrees  of  freedom  of  most  other  elements  have  transformed  as  vectors.  Cowper  et  al. 

[22]  developed  a  shell  element  for  shallow  shell  applications;  some  of  its  degrees  of  freedom 
are  similar  to  the  second  type  used  here,  but  the  transformations  are  performed  in  a  two- 
dimensional  space. 

If  (5.75)  is  expanded,  then  the  second  type  of  generalized  displacements  and  forces 
could  be  transformed  in  one  operation  in  the  same  sense  as  that  required  for  the  transforma¬ 
tion  of  the  first  type  of  displacements  and  forces.  If  we  write  the  transformation  of  the 
second  type  of  generalized  displacements  as 

{v2[  =  [t2]  |V2}  (5.88) 

then  the  examination  of  the  expansion  of  (5.75)  reveals  that  the  transformation  matrix  in 
(5.88)  will  be  in  the  form 
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(5.89) 


al  1 [  al2 

al3[a] 

a2ita]  |  a22[a] 

a23 [a] 

a31  ^  |  a32 

a33  M 

where  the  matrix  a  is  the  array  of  the  second-order  transformation  array  of  (5.66).  Now  com¬ 
bining  (5.86)  and  (5.88),  the  tr£r.sfo  nation  of  all  of  the  generalized  displacements  at  nodal 
point  i  takes  the  form 


r  t1  |  0 

n  l  *2 

L  U  i  t 


(5.90) 


Using  (5.90),  the  transformation  of  all  the  generalized  displacements,  i.e.,  degrees  of  freedom, 
of  the  element  can  be  obtained  as 


or 


_m  | 

0 

0 

Jill 

0  ! 

0 

0 

JL1 

m  | 
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Starting  with  equation  (5.88),  a  parallel  development  of  the  matrix  expressions  could 
be  given  for  the  transformation  of  the  generalized  forces.  The  final  equation,  which  has  the 
same  meaning  as  equation  (5.92),  will  simply  be  written  as 

{f[=[T]{F(  (5.93) 

where  the  transformation  matrix  is  the  same  as  in  (5.92). 

The  transformation  of  the  element  stiffness  matrix  from  the  local  element  system  to 
the  global  system  follows  directly  using  (5.92)  and  (5.93).  The  stiffness  equation  can  be 
expressed  as 

[k]M=|ff  (5-94) 

where  the  matrix  k  is  the  element  stiffness  in  local  coordinates.  Now  substituting  (5.92)  and 
(5.93)  into  (5.94),  we  obtain 

[k]  [T]  |vf  =  [TJ  |F(  (5.95) 

Then  operating  on  both  sides  of  the  above  equation  by  the  inverse  of  the  transformation 
matrix,  we  obtain  (since  orthogonal  transformation) 
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(5.96) 


[T]-1  [k]  [T]  M  =  [T]T  [k]  [T]  jvf 
The  stiffness  in  global  coordinates  can  now  be  written  as 

[K]  =  [T]T  [ k ]  [T]  (5.97) 

and  the  stiffness  equation  in  global  coordinates  can  be  expressed  as 

[K]  jv}  =  {F}  (5.98) 

The  stiffness  equation  in  local  coordinates  given  in  (5.94)  is  equivalent  to  the  incre¬ 
mental  stiffness  equation  of  (2.32).  The  notation  has  been  simplified  here  to  clarify  the 
derivation  of  the  transformations.  Thus  (5.98)  is  equivalent  to  the  incremental  equation 
(2.32)  expressed  in  global  coordinates. 

5.8  CONVERGENCE  REQUIREMENTS 

In  the  development  of  any  finite  element,  it  is  essential  that  particular  attention  be 
paid  to  satisfying  convergence  requirements.  It  is  important  that  the  convergence  be  mono¬ 
tonic  so  that  the  investigator  will  always  know  on  which  side  of  the  correct  solution  his 
current  solution  lies. 

During  the  early  application  of  the  finite  element  method,  an  unfortunate  situation  de¬ 
veloped  when  it  was  found  that  better  “convergence”  rates  were  often  obtained  if  the  known 
correct  requirements  for  convergence  were  relaxed.  These  practices  are  dangerous  and  should 
be  discouraged,  particularly  when  one  starts  to  apply  the  same  methods  to  nonlinear  problems 
where  the  solutions  are  much  more  difficult  than  in  linear  problems.  The  use  of  such 
practices  serves  only  to  cloud  the  picture. 

Irons  and  Draper  [62]  state  three  requirements  as  necessary  to  ensure  convergence: 

a.  The  assumed  displacement  functions  must  be  continuous  within  the  element  and 
the  displacements  must  be  compatible  between  adjacent  elements. 

b.  The  displacement  functions  must  include  the  rigid  body  displacements. 

c.  The  displacement  functions  must  include  the  constant  strain  states  of  the  element. 
The  second  requirement  is  actually  a  special  case  of  the  third  under  certain  conditions.  The 
first  two  of  the  above  requirements  were  mentioned  in  Chapter  1  as  two  of  the  four  primary 
difficulties  that  have  impeded  the  derivation  of  curved  shell  elements  [15]. 

The  compatibility  requirements  for  a  particular  element  come  from  the  geometrical 
boundary  conditions  obtained  in  the  derivation  when  a  variational  principle,  such  as  the 
principle  of  virtual  work  is  used.  The  necessary  boundary  conditions  for  the  shallow  shell 
theory  being  used  are  compatibility  of  the  u,  v,  and  w  displacements  and  the  derivative  of  the 
normal  displacement  w  in  the  direction  normal  to  the  edge  [38].  These  are  boundary  con¬ 
ditions  that  are  imposed  on  an  element  by  the  surrounding  elements. 
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The  continuity  of  the  u,  v,  and  w  displacements  is  ensured  since  we  are  using  equal 
order  bicubic  functions  for  each  of  the  three  displacements.  Continuity  could  not  be  guar¬ 
anteed  if  we  were  using  lesser  order  polynomials  for  u  and  v  than  for  w. 

The  continuity  of  the  rotational  displacement  is  effected  through  the  use  of  the  three 
surface  vectors.  The  normal  rotation  is  computed  through  the  use  of  equation  (5.31)  and  was 
explained  in  Section  5.3.  It  is  noted  that  the  midside  node  in  elements  1  and  2  in  Figure  5.4 
is  actually  two  different  material  points  in  the  X-Y  plane  of  each  of  the  two  elements  (not 
shown)  which  is  mapped  to  the  same  spatial  point  J  on  the  shell  surface. 

The  intent  of  the  second  requirement  for  convergence  is  to  induce  no  strain  in  the 
element  when  it  undergoes  a  rigid  body  motion.  This  requirement  on  the  assumed  displace¬ 
ments  is  strictly  in  the  context  of  the  small  displacement  theory.  Considering  equation  (5.14), 
we  see  that  all  of  the  terms  of  the  strain  tensor  must  be  included  to  describe  the  unrestricted 
motion  of  a  line  segment  undergoing  a  rigid  body  motion.  Thus  if  we  are  to  truncate  the 
expressions  for  strain  and  retain  only  the  linear  terms,  the  concept  of  inducing  strains  by  rigid 
body  motion  becomes  a  relative  matter  which  is  limited  to  the  small  displacement  theory. 

It  is  known  that  the  inclusion  of  the  rigid  body  modes  in  the  assumed  displacement 
pattern  gives  improved  linear  results  [261.  The  assumed  displacements  used  here  for  the  formu¬ 
lation,  equations  (5.24)  and  (5.25),  do  include  rigid  body  translation  and  rotation  (linear 
terms).  In  view  of  the  above  remarks  on  the  truncation  of  the  expression  for  strain,  some 
additional  information  is  needed  on  the  adequacy  of  the  formulation  for  application  to  linear 
problems. 

Examination  of  the  eigenvalues  of  the  stiffness  matrix  is  a  convenient  method  to  de¬ 
termine  the  adequacy  of  an  element  for  satisfying  the  rigid  body  requirement.  It  is  known 
that  there  will  be  zero  eigenvalues  equal  in  number  to  the  number  of  rigid  body  degrees  of 
freedom  [63]. 

During  the  generation  of  the  stiffness  matrix  for  the  quadrilateral  element,  it  was 
found  necessary  to  place  a  1.0  on  the  diagonal  for  those  degrees  of  freedom  which  are 
derivatives  with  respect  to  Z.  There  are  a  total  of  12  of  these  degrees  of  freedom.  Thus  for 
an  ascending  ordering  of  the  eigenvalues,  the  nineteenth  is  the  first  one  that  can  be  associated 
with  the  deformation  of  the  element  since  the  first  six  are  attributed  to  the  rigid  body  motion. 
We  will  subsequently  examine  the  eigenvalues  of  the  stiffness  matrices  of  several  elements  in 
order  to  demonstrate  the  above  characteristics. 
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CHAPTER  6 

IMPLEMENTATION  OF  THEORY 


The  previously  described  theory  has  been  incorporated  in  an  early  version  of  the 
general  purpose  program  developed  by  Marcal  and  his  associates  at  Brown  University.  The 
methods  of  solution  in  the  program  have  been  described  [64,  65  |  in  a  number  of  previous 
publications  and  the  basic  techniques  used  will  be  reviewed  only  briefly  here.  The  description 
will  be  primarily  limited  to  the  implementation  of  the  quadrilateral  element  and  the  previously 
described  incremental  formulation  developed  herein. 

6.1  LINEAR  ELASTIC  LARGE  DEFORMATION  SOLUTIONS 

The  study  of  the  large  deflection  of  structures  where  the  material  remains  elastic  re¬ 
quires  much  less  computational  effort  than  when  the  material  of  the  structure  is  carried  into 
the  inelastic  range.  The  [D]  matrix  of  equation  (2.30)  remains  constant  throughout  the  load¬ 
ing  history. 

Each  of  the  component  terms  in  equation  (2.30)  is  integrated  over  the  volume  and 
surface  of  the  element  in  the  local  coordinate  system.  The  region  of  integration  for  the 
element  is  shown  in  Figure  6.1.  Each  of  the  integrals  is  computed  by  the  approximate 
method  whereby  the  region  is  subdivided  into  16  subtriangles  and  the  integrand  is  then  evalu¬ 
ated  at  the  centroid  of  the  subtriangle.  To  integrate  the  function  f  over  the  element,  the 
approximation  is  made  that 


/' 


f(x,Y)  dA 


X]  f:(XC,YC)AA: 
i=l  ,16  1  1 


(6.1) 


where  the  superscript  c  indicates  the  function  evaluated  at  the  subtriangle  centroid.  Thus  the 
product  of  the  terms  appearing  in  the  integrand  of  each  of  the  respective  terms  is  evaluated  at 
the  centroid,  multiplied  by  the  contributing  area,  and  then  summed. 

The  evaluation  of  the  B1  and  B2  matrices  given  in  Table  5.1  for  the  quadrilateral  is 
required  for  the  computation  of  the  element  stiffness  matrix  and  the  equilibrium  correction 
vector  in  (2.30).  In  Table  5.1,  the  terms  in  parentheses  shown  with  bars  are  the  local  slopes 
of  the  element  surface  derived  in  Chapter  4.  The  remaining  terms  in  the  parentheses  are,  first, 
the  effect  of  initial  displacements  as  termed  by  Marcal  [6]  and,  second,  the  estimate  of  the 
displacement  that  will  occur  in  going  from  the  intermediate  configuration  to  the  current  con¬ 
figuration.  The  latter  are  prefixed  by  the  del  symbol.  The  initial  displacements  are  the  total 
displacements  that  have  occurred  up  through  the  intermediate  configuration.  The  determina¬ 
tion  of  the  estimated  displacements,  i.e.,  the  displacements  that  are  to  be  computed  for  the 
current  increment,  was  described  following  the  derivation  of  equation  (2.30)  but  additional 
detail  will  be  given  subsequently. 
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Thus,  it  is  relatively  straightforward  to  form  the  stiffness  matrix  and  the  equilibrium 
load  correction  vector  of  equation  (2.30).  The  explicit  evaluation  of  the  B1  and  B2  matrices 
at  each  of  the  integration  points  is  obtained  first  by  computing  the  derivatives  of  the  r  matrix, 
defined  in  equation  (5.25),  through  the  use  of  equation  (5.28).  The  initial  displacements  and 
the  estimated  displacements  at  the  integration  points  within  the  element  are  dependent  on  the 
initial  and  estimated  generalized  displacements  at  the  element  nodal  points.  For  the  quadri¬ 
lateral,  the  initial  and  estimated  displacements  are  in  the  w  component  only.  Their  relation¬ 
ship  to  the  generalized  displacements  at  the  nodes  is 

I^^IOiOlM^ivI  (6.2, 

a= 1.2 

which  is  written  in  the  local  element  coordinate  system.  Thus  the  initial  and  estimated  dis¬ 
placements  are  computed  for  the  integration  point  and  added  to  the  local  slope  of  the  shell  at 
the  point.  This  sum  is  then  used  in  forming  the  B1  and  B2  matrices. 

The  analysis  is  started  by  applying  a  small  portion  of  the  total  load  and  computing  the 
linear  elastic  response.  For  this  part  of  the  analysis,  the  initial  displacements  are  zero  and  the 
estimated  displacements  are  not  used  because  (a)  it  is  assumed  that  the  load  is  applied  to  the 
undeformed  structure  and  (b)  the  applied  load  is  so  small  that  the  response  is  linear.  At  this 
point  the  B1  and  B2  matrices  are  identical  and  the  computed  stiffness  matrix  is  symmetric. 

The  stiffness  matrices  are  assembled  into  the  global,  i.e.,  the  total  stiffness  matrix,  and  the 
resulting  equations  are  solved  for  the  generalized  displacements  at  the  nodes  resulting  from 
the  application  of  the  first  loading.  The  computation  of  the  consistent  loading  vectors  for 
distributed  loading  is  straightforward  and  was  described  in  Chapter  5. 

The  nonlinear  analysis  begins  at  this  point.  The  previously  computed  generalized  dis¬ 
placements  at  the  nodes  are  first  used  to  compute  the  estimated  generalized  nodal  point 
displacements  that  will  occur  in  the  next  increment  of  loading.  This  estimate  is  obtained  by 
a  simple  linear  extrapolation.  Thus  by  using  the  previous  notation  for  the  generalized  nodal 
displacements  in  the  global  system,  we  may  write  the  estimated  displacements  for  the  next 
increment  of  loading  as 


lAvwi=flAXl 

where  f  is  the  ratio  of  the  applied  loading  in  the  next  increment  to  the  applied  loading  in  the 
increment  just  computed. 

The  element  stiffness  matrices  are  now  recomputed  to  incorporate  the  nonlinear  effects. 
The  generalized  nodal  displacements  and  the  estimated  nodal  displacements  that  will  occur  in 
the  current  increment,  obtained  through  equation  (6.3),  are  first  transformed  to  the  local  ele¬ 
ment  coordinate  system  by  equation  (5.92).  The  element  internal  initial  displacements  and 
estimated  displacements  are  now  computed  for  each  integration  point  by  using  equation  (6.2). 
These  values  are  then  used  in  forming  the  B1  and  B2  matrices  required  for  the  element  stiffness 
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matrix  and  the  equilibrium  load  correction  vector  of  equation  (2.30).  The  generalized  stresses 
computed  at  the  conclusion  of  the  initial  linear  analysis  are  also  used  in  forming  the  equilib¬ 
rium  load  correction  vector. 

The  newly  computed  stiffness  matrices  are  now  transformed  and  assembled  again  into 
the  global  system.  The  total  load  vector  is  augmented  first  by  the  increment  of  direct  loading 
and  then  by  the  equilibrium  load  correction  vector  and  the  estimated  change  in  surface  loading 
due  to  a  change  in  geometry -if  this  effect  is  being  included  in  the  analysis.  These,  too,  are 
transformed  to  the  global  system.  The  estimated  displacements  within  the  element  for  use  in 
equation  (5.65)  for  computing  the  load  due  to  a  change  in  geometry  are  computed  identically 
to  those  that  were  computed  to  form  the  B1  and  B2  matrices.  The  solution  of  the  stiffness 
equations  now  yields  the  increment  of  generalized  displacements  at  the  nodes  resulting  from 
the  application  of  the  increment  of  load. 

The  incremental  displacements  of  the  nodes  that  are  computed  will  not  be  exact  since 
the  estimated  incremental  displacements  that  were  used  in  forming  the  stiffness  equations  were 
not  exact.  This  is  true  in  general  since  a  linear  extrapolation  is  used  to  obtain  the  displace¬ 
ments  and  we  know  from  the  outset  that  the  load -deflection  curve  is  nonlinear. 

An  important  advantage  of  the  method  being  described  is  that  we  know  that  the 
solution  technique  could  give  the  exact  answer  if  we  were  able  to  use  the  exact  increment  of 
displacements  in  forming  the  stiffness  equation;  this,  of  course,  is  impossible  since  it  is  what 
we  are  after  in  the  first  place.  The  word  “exact”  is  used  in  the  context  that  we  are  working 
with  the  full  equilibrium  equations  and  use  no  linearization  approximation  (except  for  estima¬ 
ting  displacement)  during  their  derivation.  Barring  divergence,  we  know  that  the  increment  of 
displacements  that  we  have  just  computed  will  give  a  better  estimate  of  the  true  incremental 
displacement  than  will  the  linearly  extrapolated  values  used  at  the  beginning  of  the  increment 
Thus  the  computed  incremental  displacements  may  now  be  used  to  form  the  revised  element 
stiffness  matrices,  the  equilibrium  load  correction  vector,  and  the  loading  vector  due  to  the 
change  in  geometry.  This  is  done  all  within  the  current  increment  of  loading,  that  is,  without 
changing  the  load.  The  process  is  repeated  until  the  desired  accuracy  is  reached. 

The  decision  to  stop  the  iteration  in  any  one  increment  of  loading  is  based  on  the 
computed  displacement  at  a  preselected  nodal  point.  If  the  newly  computed  displacement 
divided  by  the  estimated  displacement  or  the  displacement  computed  in  the  previous  iteration 
is  less  than  a  chosen  ratio,  then  the  process  is  stopped  for  the  current  increment  and  the 
analysis  goes  on  to  the  next  increment  of  loading.  Typical  values  of  the  ratio  and  the  strategy 
used  in  performing  an  analysis  will  be  described  later  when  we  consider  particular  case  studies. 

6.2  INELASTIC  LARGE  DEFORMATION  SOLUTIONS 

The  analytical  process  used  when  the  material  in  the  structures  becomes  nonlinear  and 
path  dependent  is  somewhat  similar  to  that  just  described,  but  it  is  more  involved.  We  must 
now  consider  the  changing  D  matrix  in  equation  (2.30),  which  is  written  there  for  the  general 
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continuum  case.  For  the  specialized  case  of  the  quadrilateral  shell  finite  element,  where  the 
equilibrium  equations  have  been  redefined  in  terms  of  the  generalized  stresses  and  strains 
convenient  for  plate  and  shell  problems,  the  D  matrix  takes  the  form  of  the  integral  in  equa¬ 
tion  (5.48).  At  any  point  through  the  shell,  the  stress  and  strain  must  satisfy  equation  (3.49) 
so  we  must  show  how  this  equation  is  implemented  for  the  quadrilateral  to  form  the  relation¬ 
ship  between  generalized  stress  and  strain. 

Figure  6.2  shows  a  section  through  the  thickness  of  the  element  and  the  1  1  stations 
used  to  perform  the  numerical  integration.  These  1  1  stations  are  used  at  each  of  the  1 6 
integration  points  shown  in  Figure  6.1 

Equation  (5.48)  is  rewritten  here  for  convenience 


!  _Z[p:l_~|  dZ  f  AE  ^ 
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and  the  nomenclature  is  as  defined  previously.  The  submatrix  p  is  the  stress-strain  relation¬ 
ship  given  in  parentheses  in  equation  (3.49).  For  the  present  case,  it  is  three  by  three  for 
plane  stress  in  conformance  with  the  Kirchoff  assumptions.  To  simplify  the  description,  the 
equation  will  be  rewritten  as 


(6.4) 


The  components  of  each  of  the  submatrices  are  computed  approximately  as  follows: 
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where  the  components  of  p  are  determined  at  each  step  of  the  loading.  The  integral  is  evalu¬ 
ated  by  using  the  Simpson  rule  for  approximation. 

The  total  numerical  integration  scheme  for  computing  the  element  stiffness  matrix  may 
now  be  written  as 
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where  the  symmetric  part  is  taken  as  the  stiffness  and  the  nonsymmetric  part  is  treated  as  an 
equivalent  force  as  described  in  Chapter  2. 
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The  varying  nonconservative  behavior  of  inelastic  materials  requires  recording  the  trace 
of  each  material  point  in  the  stress  and  strain  space.  Thus  for  each  of  the  1 1  stations  through 
the  thickness  at  each  of  the  16  integration  points  of  the  element,  sufficient  information  has  to 
be  recorded  to  allow  determination  of  the  current  incremental  relationship  between  stress  and 
strain.  The  three  current  stresses,  two  direct  and  one  in-plane  shear,  and  the  current  equivalent 
plastic  strain  must  be  stored  for  each  point  for  both  the  isotropic  and  kinematic  hardening 
rules.  Additionally,  the  isotropic  rule  uses  the  current  equivalent  yield  stress  whereas  the  kine¬ 
matic  rule  uses  the  center  of  the  current  yield  surface  to  determine  the  appropriate  incremental 
relationship.  The  following  description  applies  primarily  to  the  isotropic  hardening  rule. 

The  program  allows  the  option  of  scaling  the  linear  solution  to  the  load  which  causes 
yielding  at  the  highest  stressed  point  within  one  of  the  elements  or  the  option  of  incrementing 
the  load  (as  described  in  the  previous  section)  with  elastic  material  properties  before  the  plas¬ 
ticity  analysis  started.  This  second  option  is  important  for  structures  in  which  the  effects  of 
nonlinear  geometry  may  be  important,  e.g.,  shells. 

After  each  increment  of  load  is  applied  and  a  solution  is  obtained,  each  point  in  each 
element  is  tested  to  determine  whether  it  will  go  into  the  plastic  range  on  application  of  the 
next  increment  of  load.  This  is  found  from  the  estimated  generalized  nodal  point  displace¬ 
ments  described  in  the  previous  section.  Thus,  using  the  relationship  between  incremental 
generalized  strain  and  incremental  generalized  nodal  point  displacement  in  local  element 
coordinates  as 


{MHB2]tA]Me  (6.9) 

we  may  obtain  the  strain  at  any  point  through  the  thickness  by  using  equation  (5.45).  Once 
the  estimated  strain  is  computed,  the  estimated  increment  of  stress  that  will  occur  is  computed 
in  accordance  with  equation  (5.43)  where  only  elastic  material  properties  are  used  in  the 
equation.  This  estimate  of  stress  is  then  added  to  the  stress  that  exists  at  the  point.  If  on  the 
basis  of  the  von  Mises  criterion  it  is  determined  that  the  point  will  remain  elastic,  then  compu¬ 
tation  of  the  stiffness  can  continue.  If  not,  the  material  at  the  point  will  be  in  the  transition 
region  and  the  estimated  reduction  in  stiffness  will  have  to  be  determined. 

The  approach  introduced  by  Marcal  and  King  [66]  for  determining  the  incremental 
relationship  between  stress  and  strain  is  briefly  outlined  as  follows: 

a.  For  simplification  here,  it  is  assumed  that  the  increment  of  load  will  allow  the 
point  to  remain  elastic  during  part  of  the  increment,  say  the  ratio  m,  and  that  the  remaining 
part  will  be  inelastic,  (1-m).  McNamara  [67]  described  the  use  of  this  method  for  determin¬ 
ing  the  value  of  m  in  connection  with  problems  of  transient  vibration.  Figure  6.3  gives  a 
representation  of  the  method.  The  value  of  m  is  obtained  by  writing  the  von  Mises  yeild 
condition 


J2  (S°  +  m  AS)  =  1/3  0^ 


(6.10) 
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Figure  6.3  —  Calculation  of  Ratio  m  and  Normal  to  Yield  Surface 
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where  S°  is  the  stress  state  at  the  intermediate  configuration  and  AS  is  the  increment  of  stress 
that  results  if  the  increment  of  load  allows  the  material  to  remain  elastic.  A  quadratic  eaua- 
tion  in  m  results  and  the  root  is  taken  that  causes  the  stress  vector  to  be  pointed  outward 
from  the  loading  surface. 

b.  With  the  value  of  m,  the  estimated  increment  of  stress  occurring  up  to  yield  is 
added  to  the  current  stress.  The  deviatoric  part  of  these  stresses  is  then  computed. 

c.  Now  by  using  the  estimated  deviatoric  stresses,  the  total  estimated  strains  (from  the 
estimated  nodal  point  generalized  displacements  and  equation  (6.9)),  the  slope  of  the  equiva¬ 
lent  stress-equivalent  plastic  strain  curve,  and  the  equivalent  stress  in  equation  (3.50),  the 
estimated  increment  of  equivalent  plastic  strain  can  be  computed.  If  this  is  the  first  plastic 
loading  of  the  point,  the  initial  slope  of  the  equivalent  stress-equivalent  plastic  strain  curve  is 
used  in  the  computation.  The  strain  is  then  reduced  by  the  value  of  m.  Writing  the  estimated 
plastic  strain  in  finite  increment  form,  we  have 


Aef  =(l-m) 


2/3  ae  Q| j  Cjjrs 


da 

4/9  7"e  ae+-CTijCijkCffkC 
fleP 


Aers 


(6.11) 


d.  The  above -computed  equivalent  plastic  strain  allows  determination  of  the  estimated 
slope  of  the  material  curve  at  the  end  of  the  increment.  This  slope  is  used  to  determine  a 
weighted  mean  slope  that  is  assumed  to  be  acting  over  the  increment. 

e.  The  final  incremental  stress-total  strain  relationship  is  computed  to  be 


Aaij  =  1  rnCijrs —  ( 1  —  m) 


^ijkfi  gk£  gpq  Cpqrs 
da 

4^de®  °e  +aPqCP9rs<7rs 
P 


(6.12) 


f.  If  this  is  the  last  iteration  of  a  particular  load  increment,  the  equivalent  plastic  strain 
and  total  stress  are  now  updated  and  stored  by  using  equations  (6.1 1)  and  6.12),  respectively. 

If  more  iterations  are  needed  to  meet  the  desired  accuracy,  the  quantity  in  parentheses  in 
equation  (6.12),  which  corresponds  to  the  p  matrix  in  equations  (6.5)  through  (6.7),  is  used  to 
form  the  D  matrix  for  the  next  iteration. 

The  description  of  the  implementation  of  the  kinematic  hardening  rule  is  essentially 
the  same  as  given  above  except  that  the  loading  surface  translates;  this  requires  tracking  and 
storing  the  surface  origin.  Special  attention  has  to  be  taken  to  ensure  that  the  surface  trans¬ 
lates  along  the  normal  to  the  surface  for  the  plane  stress  case  under  consideration  here. 

Hibbitt  [11]  satisfied  this  requirement  by  recognizing  the  incompressibility  of  plastic  deforma¬ 
tions.  Thus  the  trace  of  the  translation  tensor  in  equation  (3.32)  must  be  equal  to  zero.  This 
is  satisfied  by  adding  to  the  two  direct  stress  translation  components  the  negative  of  their  sum. 
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CHAPTER  7 

NUMERICAL  EXAMPLES 


The  accuracy  of  the  previously  derived  shell  element  and  nonlinear  formulation  will 
now  be  examined.  We  will  first  consider  the  linear  applications  of  the  element,  then  study  a 
very  simple,  highly  nonlinear  rod  and  spring  problem,  and  finally  consider  nonlinear  applica¬ 
tions  of  the  shell  element. 

7.1  RIGID  BODY  MODES 

Our  method  for  studying  the  rigid  body  modes  of  the  quadrilateral  shell  element  was 
introduced  in  Section  5.8.  It  was  stated  there  that  the  linear  stiffness  matrix  should  have  six 
zero  eigenvalues,  corresponding  to  each  of  the  rigid  body  modes;  since  we  are  using  a  numerical 
approach,  the  six  eigenvalues  will  not  be  exactly  zero.  We  also  pointed  out  that  the  eigenvalue 
would  be  considered  in  ascending  order  and  that  the  nineteenth  would  be  the  first  one 
associated  with  the  deformation  of  the  element. 

The  results  of  the  study  of  the  eigenvalues  are  given  in  Figure  7.1.  The  stiffness  matrix 
of  an  element  that  represents  the  10-inch  sector  of  a  spherical  shell  is  analyzed.  Four  differ¬ 
ent  radii  of  the  sphere  are  considered.  The  matrix  subroutine  package  of  the  Argonne  National 
Laboratory  [69]  was  used  to  compute  the  values. 

The  results  shown  in  the  figure  clearly  demonstrate  the  characteristics  stated  in  Section 
5.8.  It  is  also  seen  that  as  the  element  becomes  more  and  more  shallow,  i.e.,  as  the  ratio  of 
element  side  length  to  radius  of  shell  decreases,  the  results  become  closer  to  what  is  expected. 
The  first  six  eigenvalues  were  very  small  compared  to  the  seventh.  Moreover,  there  was  a  big 
difference  between  the  eighteenth  and  nineteenth  eigenvalues. 

The  dependency  of  the  computed  eigenvalues  on  the  shallowness  of  the  element  indi¬ 
cates  that  for  application  to  small  displacement  problems,  the  adequacy  of  shell  fininte  ele¬ 
ments  to  satisfy  the  rigid  body  requirement  is  highly  dependent  on  the  strain -displacement 
relationships  that  are  used  in  the  element.  In  other  words,  as  the  element  becomes  flatter, 
the  shell  terms  become  less  significant  and  the  computed  values  approach  the  ideal  values. 

7.2  MEMBRANE  CYLINDER -STRESSES 

Examination  of  the  need  for  an  equilibrium  correction  term  in  the  incremental  formu¬ 
lation  (Section  2.2)  showed  that  it  arose  from  the  fact  that  the  stresses  did  not  satisfy  the 
equilibrium  equation  at  each  point  within  the  element  and  on  the  boundary.  Thus,  the 
equilibrium  equations  are  satisfied  only  in  an  integral  or  average  sense.  The  question  arises 
then  as  to  how  to  interpret  the  stresses  that  can  be  computed  at  each  one  of  the  integration 
points  shown  in  Figure  6.1.  In  other  words,  how  can  physically  meaningful  stresses  be  ob¬ 
tained?  One  way  would  be  to  make  the  elements  very  small,  smaller  than  required  for  an 
accurate  displacement  solution.  Of  course,  this  is  undesirable  because  the  additional  elements 
increase  computing  costs. 
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EIGENVALUES 
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Figure  7.1  -  Eigenvalues  of  Stiffness  Matrix  of  Quadrilateral  Element  for  Sector  of  Sphere 


When  applying  the  quadrilateral  element,  the  use  of  average  strains  to  compute  the 
stresses  has  been  found  a  very  convenient  approach.  Thus,  by  using  the  notation  of  Chapters 
5  and  6,  the  average  strain  may  be  written  as 

{£vHb2wai  ia?i  ,7 '> 

where 

IB2W  =  (7.2, 

i=l,16 

and  ASj  is  the  area  of  the  subtriangle  i  in  Figure  6.1.  Thus  equation  (7.2)  is  an  expression  for 
a  weighted  B2  matrix  with  which  the  average  strain  can  be  computed. 

The  analysis  of  the  membrane  cylinder  shown  in  Figure  7.2  demonstrates  the  argument 
for  using  the  average  strain  concept.  The  cylinder  is  loaded  externally  with  hydrostatic 
pressure.  It  is  restrained  from  moving  in  the  longitudinal  direction  and  thus  is  in  a  state  of 
plane  strain.  The  cylinder  is  idealized  with  three  elements  all  of  the  same  size. 

The  radial  displacement  computed  by  using  the  quadrilateral  element  agreed  within  2 
percent  of  the  exact  value.  However,  the  membrane  stresses  computed  at  the  16  integration 
points  ranged  from  25  percent  over  to  15  percent  under  the  exact  values.  The  average  mem¬ 
brane  stresses  computed  by  using  equation  (7.1)  were  4.5  percent  over  the  exact  values.  This 
is  much  better  than  obtainable  by  the  usual  practice  of  taking  the  stress  at  the  centroid  of  the 
element  (approximately  15  percent  under  the  correct  value). 

7.3  GRAVITY-LOADED  CYLINDER 

The  quadrilateral  element  was  evaluated  with  reference  to  the  cylindrical  shell  roof 
shown  in  Figure  7.3.  The  shell  was  loaded  by  its  own  dead  load;  it  was  supported  by 
diaphragms  (simple  support)  at  each  end  and  was  free  along  both  sides. 

The  roof  problem  has  been  studied  by  many  developers  of  finite  elements  for  shells 
and  has  become  a  bench  mark  for  comparing  the  accuracy  of  the  elements  reported.  In  a 
paper  comparing  the  finite  difference  and  finite  element  methods,  Forsberg  [70]  pointed  out 
the  difference  in  deflections  that  can  be  obtained  with  the  finite  element  method  in  applica¬ 
tion  to  this  problem.  He  reported  that  these  differences  can  be  attributed  to  the  shell 
equations  used  in  formulating  the  various  elements.  By  first  programming  the  shallow  shell 
and  then  the  deep  shell  equations  into  a  very  accurate  finite  difference  program,  he  was  able 
to  explain  the  differences  that  had  previously  been  reported  in  the  literature.  Forsberg  also 
indicated  [70]  that  the  analysis  which  uses  flat  plate  elements  converges  to  the  results  ob¬ 
tained  by  the  deep  shell  equations  implemented  in  the  finite  difference  program. 

The  convergence  of  solutions  to  the  displacement  obtained  at  the  center  of  the  free 
edge  is  plotted  in  Figure  7.4  for  several  of  the  most  recent  elements  published  and  for  the 
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Figure  7.4  -  Total  Number  of  Degrees  of  Freedom  (Unknowns) 


quadrilateral  element  reported  here.  The  two  results  obtained  by  the  finite  difference  shallow 
shell  and  deep  shell  study  of  Forsberg  [ 70 ]  are  also  shown.  The  abscissa  is  the  total  number 
of  unknowns. 

Ahmed  et  al.  [71]  extended  to  shell  problems  the  isoparametric  concept  used  for  the 
development  of  a  broad  class  of  finite  elements.  These  results  are  shown  in  Figure  7.4. 
Specializations  were  recently  added  to  make  the  original  formulation  more  effective  with 
reference  to  thin  shells  [72].  Several  schemes  were  studied,  and  two  of  the  more  rapidly 
converging  ones  are  included  in  Figure  7.4.  The  first  is  the  case  where  a  simple  2x2  Gauss 
point  integration  was  used  with  the  parabolic  element,  and  the  second  is  the  application  of  a 
cubic  element  with  modified  transverse  shears.  It  is  interesting  to  note  that  the  first  converges 
to  the  shallow  shell  results  and  the  second  to  the  deep  shell  results  obtained  by  Forsberg  [70]. 

The  result  obtained  by  the  rapidly  converging  element  of  Cowper  et  al.  [22]  is  also 
plotted  in  Figure  7.4.  This  element  is  a  conforming  type  and  uses  fifth  order  polynomials 
for  displacements,  which  explains  its  good  behavior.  Its  convergence  to  the  shallow  shell  dis¬ 
placement  of  Forsberg  [70]  is  as  expected  since  it  is  this  type  of  strain  displacement  equation 
which  is  used  in  global  coordinates  to  derive  the  element. 

Finally,  the  displacements  obtained  from  the  quadrilateral  element  reported  here  are 
given  in  Figure  7.4  for  increasing  grid  refinement.  Meshes  of  2  x  2,  3  x  3,  4  x  4,  4  x  5,  and 
4x7  were  used  in  the  study.  The  five  and  seven  subdivisions  were  in  the  longitudinal  direc¬ 
tion.  Since  the  Marguerre  shallow  shell  theory  [58]  is  applied  in  local  element  coordinates 
and  not  in  global  coordinates,  the  displacements  converge  to  the  deep  shell  displacement. 
Computing  time  for  the  3  x  3  case  on  the  CDC  6400  was  2.6  minutes,  approximately  three 
times  that  required  on  a  CDC  6600  computer. 

7.4  RING-STIFFENED  CYLINDER 

In  Section  5.3,  the  generalized  displacements  were  increased  to  12  per  node  to  include 
the  derivative  of  u,  v,  and  w  with  respect  to  z.  This  was  done  to  facilitate  the  transformation 
of  the  element  so  that  it  could  be  applied  to  branched  and  stiffened  shells. 

The  portion  of  the  ring-stiffened  cylinder  shown  in  Figure  7.5  was  analyzed  to  dem¬ 
onstrate  this  capability.  It  is  assumed  that  the  cylinder  is  infinite  in  length  with  periodic 
spacing  of  the  stiffening  rings.  Node  A  is  a  midbay  between  two  rings.  Since  the  behavior  is 
periodic,  the  idealization  shown  may  be  used  to  analyze  the  shell. 

At  points  where  the  shell  is  not  intersected  by  another  shell  or  stiffening  reinforce¬ 
ment,  such  as  nodal  point  B,  a  local  nodal  point  coordinate  system  is  established  as  described 
in  Section  4.4  and  the  degrees  of  freedom  that  are  derivatives  with  respect  to  the  shell  normal 
are  set  equal  to  zero.  This  is  generally  true  in  applying  the  element.  Since  these  degrees  of 
freedom  are  not  associated  with  the  stiffness  of  the  element,  they  are  eliminated.  This  is  not 
done  at  points  such  as  C  since  the  derivatives  with  respect  to  the  normal  on  the  shell  surface 
are  degrees  of  freedom  that  are  associated  with  the  stiffness  of  the  elements  in  the  web  of 
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RADIAL  DISPLACEMENTS 
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Figure  7.5  —  10-Element  Idealization  of  Ring- Stiffened  Cylinder 


the  stiffener.  Thus  there  is  a  maximum  of  nine  degrees  of  freedom  at  points  such  as  B,  and  a 
maximum  of  twelve  degrees  of  freedom  at  points  such  as  C. 

The  results  obtained  by  using  the  quadrilateral  element  in  the  ten -element  idealization 
are  tabulated  in  the  figure.  The  displacements  obtained  by  using  the  theory  presented  by 
Pulos  and  Salerno  [73]  are  also  shown.  The  maximum  difference  is  seen  to  be  approximately 

4.5  percent.  It  is  noted  that  the  study  could  just  as  well  have  been  done  with  five  elements 
since  there  is  no  variation  in  the  circumferential  direction. 

7.5  BAR-SPRING  PROBLEM 

The  previously  derived  incremental  formulation  is  applicable  to  all  nonlinear  static 
structural  problems  which  include  the  effects  of  large  deformations.  Figure  7.6  presents  a 
very  simple,  but  highly  nonlinear,  problem  that  is  good  for  purposes  of  evaluation  and  com¬ 
parison. 

The  bar  problem  has  been  studied  by  a  number  of  investigators,  most  recently  by 
Hibbitt  [11]  who  uses  the  example  to  demonstrate  the  error  that  occurs  in  summing  stresses 
in  an  incremental  solution  to  nonlinear  problems.  Previous  incremental  formulations  linear¬ 
ized  the  nonlinear  equation  in  order  to  obtain  a  solution;  this  caused  an  error  in  the  accumu¬ 
lated  stresses  since  the  strain,  and  hence  the  stresses,  change  quadratically  in  an  increment. 

The  results  of  the  linearized  formulation  are  plotted  in  Figure  7.6  along  with  the  results 
obtained  by  the  incremental  formulation  derived  herein.  The  finite  element  results  are  plotted 
against  the  exact  solution.  Note  that  the  formulation  derived  in  the  present  study  noticeably 
improves  the  computed  displacements. 

Since  the  formulation  that  was  derived  in  Chapter  2  is  complete,  an  iteration  process 
should  lead  to  very  accurate  answers  for  the  bar-spring  case.  This  was  found  to  be  true. 

The  results  of  the  study  are  given  in  Table  7.1.  The  linear  result  was  obtained  first 
with  a  1/1 0-pound  load.  The  displacement  was  then  computed  at  load  levels  that  were 
multiples  of  3  pounds  up  to  30  pounds.  The  tolerance  ratio  referred  to  in  the  table  is  de¬ 
scribed  in  the  last  paragraph  of  Section  6.1.  The  integers  shown  in  the  table  are  the  number 
of  iterations  required  to  go  from  one  load  level  to  the  next  in  order  to  obtain  the  indicated 
accuracy. 

The  entries  in  the  table  indicate  that  very  accurate  results  can  be  obtained  with  the 
formulation  if  one  is  willing  to  pay  the  price  for  computing.  There  may  be  cases  where 
extreme  accuracy  is  desired,  but  for  most  practical  problems,  the  high  number  of  iterations 
shown  for  the  tighter  tolerances  is  not  warranted.  The  important  point  is  that  the  accuracy  is 
there  if  it  is  needed. 
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TABLE  7.1 

BAR-SPRING  PROBLEM 
DISPLACEMENT  VERSUS  TOLERANCE  RATIO 


Load 

Exact 

Tolerance  Ratio 

lb 

Solution 

2.0 

1.75 

1.5 

1.25 

1.1 

1.05 

1.01 

1.001 

0.1 

3.0 

0.235363 

0.00625059 

0 

0.22458 

0 

0.22458 

0 

0.22458 

0 

0.22458 

1 

0.232882 

1 

0.232882 

2 

0.234794 

5 

0.235340 

6.0 

1.0000 

0 

0.593575 

0 

0.593575 

1 

0.681971 

1 

0.681971 

3 

0.794509 

5 

0.858871 

13 

0.966944 

26 

0.997661 

9.0 

1.76462 

1 

1.79608 

1 

1.79608 

1 

1.79130 

1 

1.79130 

1 

1 .76925 

1 

1 .76409 

1 

1 .76428 

1 

1 .76488 

12.0 

2.0000 

0 

1.94017 

1 

2.00431 

1 

2.00418 

2 

1.99996 

2 

1.99991 

3 

1 .99992 

3 

2.00026 

4 

2.00024 

15.0 

2.16169 

0 

2.17539 

0 

2.15485 

0 

2.15424 

1 

2.16286 

1 

2.16329 

2 

2.16168 

3 

2.16193 

4 

2.16189 

18.0 

2.28913 

0 

2.26873 

0 

2.28626 

0 

2.28646 

1 

2.29059 

1 

2.29062 

2 

2.28908 

3 

2.28936 

4 

2.28930 

21.0 

2.39610 

0 

2.40428 

0 

2.39155 

0 

2.39134 

0 

2.39144 

1 

2.39731 

2 

2.39604 

3 

2.39630 

4 

2.39624 

24.0 

2.48921 

0 

2.47746 

0 

2.48749 

0 

2.48765 

0 

2.48861 

1 

2.49023 

1 

2.49017 

3 

2.48939 

4 

2.48933 

27.0 

2.57222 

0 

2.57828 

0 

2.56939 

0 

2.56925 

0 

2.56878 

1 

2.57309 

2 

2.57210 

3 

2.57238 

4 

2.57231 

30.0 

2.64747 

0 

2.63890 

0 

2.64651 

0 

2.64663 

0 

2.64717 

1 

2.64823 

1 

2.64808 

3 

2.64761 

4 

2.64755 

7.6  SIMPLY  SUPPORTED  SQUARE  PLATE 

The  next  example  gives  a  further  comparison  between  the  linearized  formulation  and 
the  formulation  developed  in  this  study.  The  structure  is  a  simply  supported  square  plate 
(Figure  7.7)  subjected  to  uniform  pressure  loading.  The  plate  is  restrained  against  in-plane 
motion  and  thus  causes  the  buildup  of  membrane  forces.  The  thickness  is  0.1  in.,  Young’s 
modulus  is  30  million  psi,  and  Poisson’s  ratio  is’ 0.3 16. 

Levy  [74]  solved  this  problem  by  a  trigonometric  series  expansion.  Bergan  [47]  used 
16  elements  for  a  finite  element  solution  to  the  problem  and  performed  a  Newton-type 
iteration  at  each  load  level  to  improve  the  results.  Their  results  are  shown  in  Figure  7.7  along 
with  the  Hibbitt  results  [11]  for  the  linearized  formulation  and  the  formulation  developed  in 
the  present  study.  Both  runs  with  each  of  Hibbitt  and  present  formulations  were  performed 
with  the  tolerance  ratio  set  equal  to  2.0  in  order  to  compare  running  times.  Since  the  formu¬ 
lation  presented  in  Chapter  2  does  not  require  the  computation  of  the  initial  stress,  i.e., 
geometric  stiffness  matrix,  one  would  expect  the  running  time  to  be  less  for  this  theory.  The 
formulation  of  this  study  took  57  percent  of  the  time  required  for  the  theory  presented  by 
Hibbitt  [11].  One  would  not  expect  this  same  reduction  in  all  cases  since  as  size  the  problem 
gets  larger,  the  computation  of  stiffness  matrices  becomes  a  smaller  part  of  the  total  computing 
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O  LINEARIZED  FORMULATION 
BERGAN  (47)  16-ELEMENT 


MIDPOINT  DEFLECTION  (IN) 

Figure  7.7  —  Load  Deflection  Curve  for  Square  Plate 


time.  Note  that  again  that  the  present  formulation  gives  improved  results  over  the  linearized 
formulation  for  the  same  idealization  and  iterative  scheme. 

7.7  CENTRALLY  LOADED  SPHERICAL  CAP 

The  centrally  loaded  aluminum  cap  shown  in  Figure  7.8  was  analyzed  to  give  an 
example  of  inelastic  material  behavior.  A  90-deg  segment  of  the  cap  was  idealized  with  seven 
elements.  This  case  and  the  experimental  results  are  from  Penning  and  Thurston  [75].  The 
shell  was  fixed  at  the  boundary,  Young’s  modulus  was  equal  to  1  x  107,  and  Poisson’s  ratio 
was  0.33.  The  proportional  limit  of  the  effective  stress-effective  strain  curve  was  37,000  psi. 
To  model  the  strain -hardening  portion  of  the  curve,  three  linear  segments  with  the  following 
properties  were  used:  0.5  x  107,  0.0;  0.15  x  1 07 ,  0.9  x  10'^;  and  0.44  x  106,  0.28  x  10'^. 
The  first  number  in  each  pair  gives  the  slope  and  the  second  number  the  starting  equivalent 
plastic  strain  of  the  segment. 

The  comparison  between  experimental  results  and  the  computed  results  obtained  here 
is  shown  in  Figure  7.9.  Two  different  loading  cases  are  studied,  concentrated  load  and  the 
case  where  the  load  is  distributed  over  an  area  of  0.5 -in.  radius  at  the  apex.  The  experimental 
result  for  the  first  case  is  for  a  loaded  area  with  a  radius  equal  to  0.0625.  It  was  felt  that 
having  the  load  concentrated  at  a  point  would  be  sufficient  to  model  the  physical  behavior. 

The  considerable  difference  in  the  displacements  for  the  two  cases  shown  in  Figure  7.9 
can  be  explained  by  the  early  yielding  that  occurred  in  the  point-loaded  case.  Only  this  case 
was  run  by  both  the  linearized  formulation  and  the  formulation  presented  in  this  study.  The 
tolerance  ratio  used  in  all  runs  was  1 .2.  The  present  formulation  required  63  percent  of  the 
computing  time  required  for  the  linearized  formulation. 

7.8  PEAR-SHAPED  CYLINDER 

The  next  problem,  shown  in  Figure  7.10,  poses  a  good  exercise  for  both  the  quadri¬ 
lateral  element  and  the  nonlinear  incremental  formulation  presented  herein.  This  example  is 
taken  from  Hartung  [76]  whose  results  were  obtained  with  the  finite  difference  method.  The 
structure  is  similar  to  a  proposed  configuration  for  a  space  shuttle  fuselage. 

The  load  was  applied  to  the  end  of  the  cylinder  by  a  uniform  end  shortening.  The 
cylinder  was  modeled  by  eight  of  the  quadrilateral  shell  elements.  Only  one-fourth  of  the 
total  area  of  the  cylinder  had  to  be  modeled  because  of  the  double  symmetry.  The  tolerance 
ratio  used  in  the  analysis  was  1.05.  Again,  this  is  the  ratio  described  in  the  last  paragraph  of 
Section  6.1.  The  integers  in  parentheses  next  to  the  data  points  indicate  the  number  of 
iterations  required  to  go  from  one  load  step  to  the  next.  Thus  in  the  first  increment  of  load¬ 
ing  after  the  linear  solution,  i.e.,  in  going  from  100  to  200  lb  of  loading,  the  equations  were 
formed  the  one  required  time  plus  five  additional  times.  In  the  upper  limits  of  the  loading, 
the  equations  were  merely  assembled  and  solved  one  additional  time. 
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Figure  7.8  —  Centrally  Loaded  Spherical  Cap 
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Figure  7.9  -  Comparison  of  Computed  and  Experimental  Results  for 
Centrally  Loaded  Spherical  Cap 
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Figure  7.10  —  Load-Displacement  for  Pear-Shaped  Cylinder 
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The  Hartung  analysis  [76]  predicted  a  collapse  at  about  2400  lb  whereas  the  present 
analysis  was  able  to  determine  an  equilibrium  position  at  2800  lb,  the  load  at  which  the 
computation  was  stopped.  The  results  given  here  for  the  finite  element  analysis  cannot  be 
regarded  as  conclusive  until  the  same  result  is  obtained  in  another  analysis  performed  with  a 
greater  number  of  elements.  This  additional  run  would  have  to  be  made  if  the  present 
analysis  were  being  done  for  production  or  design  work.  At  this  point  we  note  only  that  it 
is  very  encouraging  that  a  reasonable  analysis  is  obtainable  with  only  eight  elements  and  that 
the  maximum  load  is  reached  in  eight  increments  of  loading. 
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CHAPTER  8 

SUMMARY  AND  RECOMMENDATIONS  FOR  FURTHER  WORK 


8.1  SUMMARY 

The  first  part  of  this  dissertation  presented  a  derivation  for  an  incremental  solution 
procedure  for  general  structural  mechanics  problems.  It  was  shown  that  the  contribution  of 
the  sum  of  all  previous  increments  had  to  be  included  in  the  solution  for  each  load  increment. 
The  resulting  nonlinear  equilibrium  equations  were  then  manipulated  and  put  into  a  form 
which  allowed  a  combined  incremental-iterative  approach  to  their  solution.  The  finite  ele¬ 
ment  formulation  was  then  used  to  put  the  equilibrium  equations  into  a  form  for  their 
numerical  solution.  Since  the  resulting  set  of  equations  in  the  displacements  turned  out  to  be 
nonlinear  and  nonsymmetric  in  the  coefficients,  a  solution  procedure  was  introduced  which 
uses  present  linear-symmetric  equation  solvers.  The  result  of  all  the  above  was  a  direct  solu¬ 
tion  procedure  whereby  the  nonlinear  incremental  equations  can  be  solved  by  using  existing 
solution  procedures  without  the  need  to  first  linearize  the  equations. 

Next,  the  basic  theory  of  plasticity  required  for  the  stress-strain  relationships  was 
reviewed.  The  description  considered  only  small  strain  problems  which  may  be  assumed  to  be 
independent  of  temperature  and  time.  Although  this  was  not  a  primary  part  of  this  disserta¬ 
tion,  its  inclusion  demonstrated  the  versatility  of  the  previously  derived  formulation  of  the 
incremental  equilibrium  equations  and  solution  procedure.  The  inclusion  of  the  plastic  ma¬ 
terial  behavior  was  possible  because  a  linear  stress  relationship  was  obtained  between  an  incre¬ 
ment  of  total  stress  in  terms  of  an  increment  of  total  strain,  and  this  is  what  was  needed  in 
the  derivation  of  the  incremental  formulation. 

The  geometry  of  a  general  quadratic  surface  was  considered  next  in  preparation  for  the 
derivation  of  the  quadrilateral  shell  element.  The  development  work  done  here  provides  the 
necessary  tools  for  performing  an  analysis  of  a  structure  which  may  be  made  up  of  several 
different  quadratic  surfaces.  The  slopes  of  the  surface  in  terms  of  the  local  system  can  be 
computed  by  using  a  coordinate  transformation  in  which  the  quadratic  form  of  the  surface  is 
expressed  in  terms  of  the  variables  of  each  individual  element  coordinate  system.  This  then 
enables  the  exact  shell  geometry  to  be  used  for  computing  the  element  stiffness  matrices  and 
for  correcting  the  load  to  the  base  plane  projected  reference  system.  The  work  on  geometry 
was  concluded  with  emphasis  on  the  surface  vectors  which  were  to  be  used  later  in  satisfying 
the  element  convergence  requirements  and  in  performing  transformations. 

The  derivation  of  the  quadrilateral  shell  element  was  performed  using  the  Marguerre 
shallow  shell  strain-displacement  relationships  [58],  The  assumed  displacement  distribution 
used  in  the  element  is  the  same  as  that  first  proposed  by  Fraeijs  de  Veubeke  [59]  except  that  a 
linearization  of  the  derivative  degrees  of  freedom  was  imposed  along  the  boundary  to  eliminate 
the  midside  nodes.  The  previous  incremental  formulation  of  the  equilibrium  equations  was 
then  specialized  for  the  quadrilateral  element  in  terms  of  the  generalized  stresses  and  strains; 
this  resulted  in  expressions  for  the  matrices  used  in  computing  the  stiffness  matrix  and  the 
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equilibrium  correction  vector.  The  form  of  the  load  vectors  due  to  pressure  loading,  the 
change  in  loading  due  to  a  change  in  geometry,  and  the  body  force  loading  were  derived. 
Transformations  required  for  the  element  stiffness  matrix  and  its  degrees  of  freedom  were  then 
developed.  These  consistent  transformations  eventually  facilitated  the  study  of  branched  and 
stiffened  shells.  The  derivation  of  the  quadrilateral  element  was  concluded  with  an  explana¬ 
tion  of  the  work  done  to  ensure  that  the  element  satisfied  the  necessary  convergence  require¬ 
ments. 

The  vehicle  used  for  implementing  all  of  the  above  and  without  which  the  work  would 
have  been  practically  impossible  is  the  general  purpose  nonlinear  finite  element  computer  pro¬ 
gram  developed  by  Marcal  and  his  associates  at  Brown  University.  The  format  of  the  program 
is  such  that  it  is  a  very  valuable  research  tool,  and  it  is  this  capability  that  was  exploited  in 
performing  the  present  work. 

The  primary  work  on  the  dissertation  was  concluded  with  the  application  of  the 
quadrilateral  element  and  the  nonlinear  incremental  formulation.  Examination  of  the  adequacy 
of  the  shell  element  to  satisfy  the  rigid  body  displacement  requirement  indicated  that  as  the 
element  became  shallower,  it  had  an  increased  ability  to  undergo  rigid  body  motion  without 
induced  strains.  The  examination  of  the  stresses  in  the  element  strengthened  past  experiences 
with  higher  order  finite  elements  which  had  indicated  that  a  weighted  mean  or  average  of  the 
stresses  provides  far  better  values  that  have  more  physical  meaning.  The  linear  application  of 
the  shell  element  appears  to  work  well  for  the  cases  studied,  e.g.,  in  application  to  a  gravity  - 
loaded  cylinder,  it  gave  very  good  results  compared  with  previously  published  elements.  Its 
application  to  a  simple  bar-spring  problem  proved  to  be  very  interesting  and  on  the  basis  of 
agreement  achieved  with  exact  values,  the  formulation  appears  to  have  the  capability 
to  produce  very  accurate  results.  The  use  of  the  quadrilateral  element  in  conjunction  with  the 
incremental  formulation  introduced  here  gives  improved  behavior  over  the  previous  linearized 
formulation.  The  application  of  the  element  and  the  formulation  to  a  pear-shaped  cylinder 
provided  very  encouraging  results.  This  is  the  only  known  finite  element  solution  to  this 
problem. 

8.2  RECOMMENDATIONS  FOR  FURTHER  WORK 

There  are  a  number  of  possible  extensions  to  the  work  reported  in  this  dissertation. 

The  most  rewarding  would  be  efforts  to  reduce  the  inordinate  amount  of  computing  time  re¬ 
quired  for  the  nonlinear  analysis  of  structures. 

The  nonlinear  incremental  formulation  introduced  herein  uses  a  linear  extrapolation 
method  for  estimating  the  amount  of  displacement  that  will  occur  in  going  from  the  inter¬ 
mediate  configuration  to  the  current  configuration.  Some  higher  order  scheme  should  be 
investigated  to  determine  whether  it  is  possible  to  reduce  the  number  of  iterations  required  to 
reach  convergence. 
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Another  variation  on  the  present  incremental  formulation  that  warrants  investigation  is 
the  treatment  of  all  nonlinearities  as  pseudo  forces.  As  noted  in  the  introduction,  the  non¬ 
linearity  due  to  inelastic  material  behavior  is  treated  this  way  in  the  initial  strain  approach. 
There  is  a  particular  advantage  in  doing  this,  particularly  for  large  problems,  since  the  system 
of  stiffness  equations  has  to  be  solved  only  once  during  the  solution  process.  (The  evaluation 
of  earlier  formulations  indicated  that  the  equations  had  to  be  solved  for  each  iteration.)  Thus, 
if  a  method  could  be  perfected  to  treat  all  of  the  nonlinearities  as  pseudo  forces,  the  equiva¬ 
lent  forces  could  be  computed  and  then  back-substituted  to  obtain  the  displacements.  How¬ 
ever,  efforts  to  treat  geometric  nonlinearities  as  pseudo  forces  have  not  been  too  successful. 

A  pilot  study  was  conducted  to  determine  whether  the  formulation  developed  herein 
was  sufficiently  accurate  and  stable  to  permit  all  of  the  geometric  nonlinearities  to  be  trans¬ 
posed  the  right  side  and  to  treat  them  as  equivalent  forces.  Using  similar  symbols  as  given  in 
equation  (2.32),  we  may  write  the  resulting  equation  as 

IK]*-  |av|  =  |p}  +  |ap|  +  |f|  +  |af[  +  |agf[  -  |e[  -  [K]N  |AVest}  (8.1 ) 

where  now  the  superscripts  on  the  K  matrix  indicate  the  linear  and  nonlinear  parts. 

Equation  (8.1)  was  applied  first  to  the  bar-spring  problem  studied  previously  (Section 
7.5).  A  tolerance  ratio  of  1.001  was  used  in  the  computation.  The  results  obtained  by  using 
equation  (8.1)  are  compared  in  Figure  8.1  with  experimental  values  and  with  the  results  of  a 
linearized  formulation  by  Haisler  et  al.[77]  for  load  increments  of  0.2  and  0.05  pounds.  The 
results  agreed  within  two  or  three  significant  digits  for  the  similar  case  given  in  Table  7.1  up 
to  a  certain  level. 

Loads  were  applied  at  the  same  levels  in  the  bar-spring  problem  as  given  in  Table  7.1, 
but  the  computation  would  not  settle  down  for  the  24-lb  load  level  and  simply  oscillated 
about  the  correct  displacement.  It  is  believed,  however,  that  this  problem  could  be  overcome 
by  improving  the  controls  used  for  determining  convergence. 

Equation  (8.1)  was  also  utilized  to  analyze  the  pear-shaped  cylinder  studied  previously 
(Section  7.8).  The  tolerance  ratio  of  1.05  was  again  used.  The  solution  process  of  equation 
(8.1)  behaved  very  well  and  results  obtained  were  identical  to  those  given  in  Figure  7.10,  in¬ 
cluding  the  number  of  iterations  to  reach  convergence. 

Finally,  it  is  recommended  that  the  incremental  formulation  presented  in  this  disserta¬ 
tion  be  extended  to  dynamic  problems  by  including  inertia  effects.  This  area  of  research  may 
prove  fruitful  inasmuch  as  the  solution  of  nonlinear  dynamic  problems  is  extremely  difficult. 
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Figure  8.1  -  Nonlinear  Geometric  Pseudo-Force  Solution 
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